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ABSTRACT 

This thesis is a contribution to the developing art of 
computer analysis for management guidance and computer control 
for economic operation of large electric generating systems. 
Particular attention is given to systems having large amounts 
of hydro storage where plant efficiencies depend on the cumula¬ 
tive effects of past storage drawdown and where present operation 
must take into account predictions of future water resources and 
system load. For this type of system, the problem is mathemat¬ 
ically formulated and procedures are developed for successively 
improving a proposed mode of operation so as to reduce an ef¬ 
fective cost of system operation, using specific predictions of 
stream flows and system load and taking into account various 
operating limitations. The semi-automatic implementation of 
these procedures is demonstrated using the Whirlwind I digital 
computer. 
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penalty cost dependent on turbine discharge at plant i 
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power loss in transmission (mw). 

total plant discharge at plant 1 (ksf). 

D^t) +J ± ( t). 
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storage elevation at dam i (ft). 

y imin = minimum allowable storage elevation at 
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tailwater elevation at dam i (ft). 
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factor determining size of correction when proceeding 
along path of steepest descent. 

incremental cost of replacing hydro deficiency 
(|/mw hr) 

modified incremental cost of replacing hydro deficiency 
($/mwhr) to take into account transmission loss. 

finite change made in 8^ at time t (ft). 

simplified plant efficiency which is constant at 
a given head (kw/efs). 

conversion factor times generator and transformer 
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Chapter I 

THE ECONOMIC COORDINATION OF POWER SYSTEMS 
1.1 Introduction 

The electric power systems are today the custodian and 
dispenser of a large portion of the nation’s energy resources. 
The economic conversion of this energy from potential to usable 
form and its economic distribution contribute to the well-being 
of the nation. The electric industry is currently witnessing 
the physical and operational integration of formerly isolated 
power systems into more effective organizations for the coordi¬ 
nation of resource-consumption over large areas (Sporn, 51) \ 
This coordination seeks to expend these resources prudently so 
as to give maximum return, or to require minimum expenditure 
for the fulfillment of the needs of the area. This requires 
the maximum exploitation of the most economical sources of 
energy, the minimization of transportation loss, the economic 
and social balance among the several benefits in multi-purpose 
river developments, and the use of minimum amounts of equip¬ 
ment, supplies and personnel. 

This marshalling of the energy resources, generating 
plants and power distribution systems of a large area into an 
integrated power system promotes economy in many ways. It 
makes possible the use of the largest units or plants justified 

"^Refers to item 51 in the bibliography, by P. Sporn. This 
notation is used throughout this treatise. 
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by the requirements of the system, regardless of the require¬ 
ments of the local area, allowing for the development of all 
resources to their maximum. It also makes possible the inte¬ 
gration of loads of diverse characteristics, thereby improving 
load factors. It makes available multi-source supplies to all 
points of the system, allowing thereby a reduction of safety 
factors in design and construction and reduction in the amount 
of capacity designated for spinning reserves. 

Thus, as the integration grows the opportunities for 
sizable operational economies and resource conservation in¬ 
crease. At the same time, however, the number of variables 
and their interrelationships also multiply so that the tasks 
of evaluating and coordinating system operations become in¬ 
creasingly complex. The precise coordination and evaluation 
of the large integrated power system become matters of great 
importance, of increasing complexity, and worthy of consider¬ 
able effort and investment. 

1.2 Types of Power System Economy Control 

First requirements of precise control are the mathematical 
formulations of the relationships among system variables, and 
the criterion for best operation of a given system with speci¬ 
fied plants, distribution system, loads and resources. The 
existence of procedures for determining optimum operation of 
a given system then makes possible the measure of best per¬ 
formance of different systems and hence, makes possible an 
evaluation of alternate proposals. 

Different types of plants and resources, in different 
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systems, require different mathematical models, somewhat 
modified criteria for best performance, and different proce¬ 
dures for determining optimum operation. For discussion 
purposes, we establish four categories of systems, which 
we call urban thermal, interconnected thermal, short-range 
hydro-thermal, and long-range hydro-thermal. 

In the first type, thermal energy is used exclusively 
and all generating plants are close to a load center. There¬ 
fore, transmission loss can be neglected and reduction of fuel 
costs is the primary consideration. The problem is the allo¬ 
cation of load among the generation plants in accordance with 
their capacities and efficiencies, so as to minimize the com¬ 
bined fuel costs with the sum of the power generations equal 
to the system load. 

Systems in the second category are characterized by 
appreciable transmission loss in long-interconnections between 
plants using exclusively thermal energy. Hence, transmission 
loss as a function of load allocation must be added to the 
mathematical model. The problem then is the allocation of 
load among the generation plants so as to minimize the com¬ 
bined fuel costs with the sum of the power generation equal to 
the system load plus the transmission loss. 

Systems employing hydro-sources in conjunction with 
thermal sources are sub-divided into two categories, long- 
range hydro-thermal and short-range hydro-thermal systems, 
both of which may or may not require the consideration of 
transmission loss. In the short-range problem, load alloca¬ 
tions are constrained so as to consume only a specified amount 
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of water at each hydro-plant over a short future time inter¬ 
val, such as one week. Variations in elevations and plant 
efficiencies during this short interval have small net effect 
on the optimum operation, so generations are still the basic 
plant variables. In systems with run of river and pondage 
plants, the specifications on the weekly amounts of water to 
be used at each plant are the anticipated availabilities 
based on short-range predictions of stream flows. In systems 
with large storage, where appreciable cumulative changes in 
plant elevations and efficiencies occur, these weekly speci¬ 
fications must be the results of a long-range optimization. 

The short range problem is the allocation of load among the 
hydro- and thermal-plants, within their capacities, so as to 
minimize the total fuel cost over a short future time interval 
with a specified average power from each hydro-plant and with 
the sum of hydro and thermal powers equal to the system loads 
plus transmission losses. 

Systems in the long-range category are characterized by 
appreciable influence of current operation on long-range 
economy due to cumulative changes in storage elevation and 
plant efficiencies. Thus in these cases operation is not 
adequately described by power generation alone, but rather by 
the particular combinations of elevation and rate of change 
of elevation which, together with natural stream flows, deter¬ 
mine generation. It follows that long-range prediction of 
stream flows and system load enter the problem explicitly. 
Furthermore, unilateral constraints imposed by various hydro 
operating limitations must be accounted for in the mathematical 
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model. In a more general case, the thermal sources are sup¬ 
plemented by imports from outside the system under study or 
contracted-for industrial curtailments in times of power 
deficiency. All such sources of energy are high-cost re¬ 
placement of hydro deficiency. The problem is the allocation 
of load among the hydro plants within their operating limita¬ 
tions so as to minimize over a long future time interval the 
cost of replacement of hydro deficiency, with proper hydraulic 
dependencies among hydro plants on interconnected streams, and 
with the sum of the power generations equal to system load 
plus transmission losses. This treatise is concerned primarily 
with systems in the long range category. 



Chapter II 

THE PRESENT STATE OF THE ART 
2.1 Introduction 

In the following a brief resume is presented of current 
practices and recent developments in economic loading of 
power systems. The past fifteen years have seen wide-spread 
interest and activity in this field. As a result system 
operation is rapidly progressing from an art toward a science. 

System economy is determined by decisions made at several 
levels of supervision. Adequate treatment of the system con¬ 
trol problem requires an understanding of the elements in¬ 
volved in economy loading at each level, the development of 
methods peculiar to the needs of problems in each level, and 
the coordination of operations at the different levels. At 
progressively higher levels of supervision in an integrated 
system are: the allocation of load among different units 
within a thermal plant; the allocation of load among thermal 
plants; the allocation of load on a short range basis, among 
hydro and thermal plants; and lastly, the timing of the use 
of water resources in hydro storage plants and the average 
load allocation between hydro and thermal sources, on a long- 
range basis. Particular progress in the study of the lower 
levels of this control hierarchy are discussed below. 
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2.2 Loading Dissimilar Units 

Various methods have been devised for allocating load 
among turbine-generator units with different heat rates 
(btu per kilowatt hour) in the same thermal plant (Steinberg 
and Smith, 56, p. $9)* Among them are (1) base loading to 
capacity, in which machines are loaded to capacity in the 
order of their over-all efficiencies, (2) base loading to 
most efficient load, which is similar except that all ma¬ 
chines are loaded in sequence according to their heat rates 
up to their most efficient loads and then to capacity in the 
same order, (3) proportionate to capacity loading, (4) loading 
proportionate to most efficient load, in which proportionate 
loading is used up to the most efficient loads of each machine 
after which loads are allocated in proportion to the differ¬ 
ence between capacity and most efficient load, (5) incremental 
loading in which all units operate at equal incremental heat 
rates. The differences in economy obtained by the various 
methods depend upon the characteristics of the units involved, 
but it has been shown that the latter, incremental loading, in 
all cases gives best economy. 

2.3 Incremental Loading 

An illustration of incremental loading, given by Steinberg 
and Smith ( 56 , p. 9), for an idealized case, is shown in figure 
2.1. The combined curve of figure 2.1b is easily constructed 

as follows: At an incremental rate of 2.0, machine A delivers 
5.0 units output and machine B delivers 2.5 units output, pro¬ 
viding a point, at an output of 7.5 and a rate of 2.0, on the 
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combined curve. Other points are similarly obtained. From 
the combined curve, then, for a given total load, the machine 
loadings giving equal incremental rates can be read off 
directly. In practice, the loading sequence is conveniently 
set up by tabulating all units at all valve settings, in the 
order of the average incremental heat rate per valve setting 
(Steinberg and Smith, 56 , p. 97; Hahn, 24-). 

In some cases, with machines of unusual characteristics, 
care must be exercised in applying the incremental loading 
criterion. In general, conditions for its applicability are 
(56, p. 7): 




(bj 


Fig. 2.1 a)Unit input-output curves. b)Incremental Rate Curves 
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a) the unit input-output curves must be continuous, 

b) the input-output curves must be well-behaved, so 
that only one cost minimum exists in the region 

of operation. This is the case if the incremental 
rate curves are monotonic increasing, even though 
they may have points of discontinuities. 

Use of the incremental heat rate takes into account unit 
efficiency and consequently fuel costs. Production costs 
other than fuel costs, particularly maintenance costs, may 
also be included as a percentage of incremental fuel costs 
based on the long-term ratio between maintenance costs and 
fuel costs. However, this refinement ordinarily does not 
appreciably influence the plant operation. 

The allocation of load among plants in an urban thermal 
system likewise is determined by incremental loading. The 
actual incremental heat rate curves for each plant have dis¬ 
continuities at the points where units are added. However, 
good solutions are obtained by using smoothed heat rate curves, 
and the corresponding smooth incremental heat rate curves 
(Steinberg and Smith, 56 ). 

Incremental loading, in effect, determines the minimum of 
a multi-dimensional surface. It does not determine what ele¬ 
ments go in to form that surface; that is, it does not deter¬ 
mine which units should be on the line for a given load. It 
determines only how to load a given set of units. The deter- 

$ 

mination of which units to be placed on the line for a given 
total plant load depends on the actual heat rates of the units 
rather than their incremental heat rates. The less efficient 
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unit is added to the line when the combined heat rates in¬ 
cluding the new unit becomes less than the combined heat 
rates without the new unit. Spinning reserv e , of course, 
modifies the number of units on the line with consequent 
sacrifice of efficiency. 

2.4 Transmission Loss Formulas 

Transmission losses are generally considered to be of 
secondary importance compared to the fuel costs or the avail¬ 
ability of water resources. For this reason, it has been 
felt that transmission losses can be taken into account 
through the use of somewhat approximate loss formulas (George, 
23; Ward, 64; Ward, Eaton and Hale, 65; Kirchmeyer and 
Stagg, 28), A compact loss formula has been developed based 
on assumptions that are not materially violated in the or¬ 
dinary operation of many systems. The basic assumptions made 
are (Kirchmeyer and Stagg, 28): 

1. The equivalent load current at any bus remains a 
constant complex faction of the total equivalent 
load current. The equivalent load current at a 
bus is defined to be the sum of the line charging, 
synchronous condenser, and load currents at the 
bus. 

2. The generator bus voltage magnitudes are assumed to 
remain constant. 

3. The ratio of reactive power to real power of any 
source is assumed to remain at a fixed value. 

The generator bus angles are assumed to remain 


4. 


constant 


The resulting transmission loss formula is: 

fl= (2-1) 

where Sm = real power output of generator M 

B nm * system constants determined usually from network 

analyzer tests 

\ 

2.5 Coordination of Fuel Costs and Transmission Loss 

To coordinate fuel costs with transmission losses one 
seeks that load allocation which will minimize total system 
fuel costs subject to the constraint that the summation of 
the plant generations shall equal the system load plus the 
transmission loss. The latter, of course, is a function of 
the load allocation. A convenient form for the necessary and 
sufficient conditions for such a minimum (or a maximum or 
saddle-point) to exist is (Kirchmeyer and Stagg, 29): 

+ A ii = A <2 - 2) 

n — 1, 2, 3, .... plants 

where F n » the fuel costs of plant n ($/hr) 

g n = the power output of plant n (mw) 

p L » the transmission loss (mw) 

A = the Lagrangian multiplier 

Often smoothed plant heat-rate curves and the corresponding 
incremental heat rates are used, neglecting discontinuities 
where additional units are placed on the line. It is then 
possible to obtain approximate analytical expressions for the 





12 


fuel costs and incremental fuel costs as follows: 


F n a f n0 * f nl S n t f n2 S n ‘ 

ila * f ♦ r g 

5 g n nl n2 B n 


($/hr) 


(|/mwhr) 


(2-3) 


(2-4) 


Also, using the above mentioned approximate expressions for 
transmission loss, the incremental transmission loss is given 


= Z 2 g i 

n i 


(2-5) 


Equation (2-2), which must be satisfied at each plant, becomes 


therefore: 

*n2 ^ £-^ni * A- F n l* n * 2 > 3,... 

( 2 - 6 ) 

Some further approximations to these conditional equations have 

also been proposed, facilitating the solution of these equations 
on the network analyzer (George, Page, and Ward, 22) or on 
standard digital computers (Kirchmeyer and McDaniels, 30). The 
errors due to these approximations were evaluated by Kirchmeyer 
and Stagg (29), and were found to be small. 


2.6 Hydro Plant Characteristics 

Typical plant efficiencies for hydro plants are shown in 
Figs. 2.2 and 2.3 (Strowger, 58, 59). In practice, loads are 
assigned so as to keep operation as much as possible near the 
points of maximum efficiency on the curves for a given head. 
The envelope of these curves may either increase or decrease 
with increased plant output. For the low-head plants improve- 



FIGURE 2.2 PLANT EFFICIENCIES FOR A LOW HEAD PLANT 
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FIGURE 2.3 PLANT EFFICIENCIES FOR A HIGH HEAD PLANT 
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ment in efficiency is obtained as units are added because the 
leakage loss in the standby units is either reduced or elim¬ 
inated. On the other hand, long pipe lines in high-head plants 
have conduit loss which varies as the square of velocity. Thus, 
as units are added in high-head plants, leakage loss is reduced, 
but conduit loss increases, producing an eventual decrease in 
efficiency. 

Ordinarily all units in a hydro-plant are identical so 
that the load is divided equally among generating units. If 
the generating units are not identical, station load should be 
divided among the units operating so as to obtain equal values 
of incremental efficiency. That station load should be carried 
which corresponds to the highest efficiency point of the station 
or at least the highest efficiency point for a given number of 
units on the line (Reid, 46 ; Schamberger, 49; Strowger, 5&, 
59). 

2.7 Hydro Thermal System Operation 

Little has been published on the economy loading problem 
involving hydro-plants (Justin and Mervine, 27; Lane, 35). 
However, certain operating practices are well established, par¬ 
ticularly in the case of run of the river and pondage plants 
concerned with short term optimization. 

The short-range hydro load allocation is constrained so 
that a specified amount of water is used at each hydro plant 
during a short future time interval, such as one week. Some 
features of the load allocation among thermal plants, run of 
the river, and pondage plants are as follows: 
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1) Gauge readings of pondage levels and stream flows are 
made periodically (perhaps twice daily) 

2) From current records of the discharges at each dam 
and the time of flow between dams, run-off between 
dams can be evaluated. 

3) River flow and run-off can then be extrapolated sev¬ 
eral days in advance with reasonable accuracy. 

(Flash floods are extrapolated differently from stream 
flow, depending upon precipitation run-off relations 
for the area.) 

4) The stream flows must be used immediately for genera¬ 
tion in run of the river plants. 

5) The total water in-flow into each pond minus the de¬ 
cided upon change in pondage level for the coming 
week determine the amount of water to be used at 

each pondage plant for the coming week. 

6) Pondage water use is distributed throughout the week 
so as to ’’fill-in” the upper portion of the week’s 
hourly load pattern. 

7) There should remain a nearly flat hydro deficiency 
which must be taken care of by thermal sources allo¬ 
cated according to incremental rates. Flat steam 

is a sign of good operating economy. 

Only an inadequate amount has as yet been published on the 
mathematical formulation of the load allocation among hydro 
plants and between hydro plants and thermal plants (Masse, 38). 
An interesting first step in this direction was taken by 
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Chandler and Gabriel in 1950 (11), in their coordination of 
water utilization with transmission loss. As an extension of 
the work of George Page and Ward (22), the problem was formu¬ 
lated on the basis of a given head at each hydro plant. In¬ 
stead of fuel consumption and transmission loss, transmission 
loss and the reduction of storage heads were minimized at any 
one time. However, it did not include the essential ingredients 
of either the short-range or long-range load allocation prob¬ 
lems, namely: a) the constraint to utilize to best advantage 
a predetermined amount of water at each plant during a coming 
short-time interval in a short-range problem, and b) the taking 
into account of the effects of cumulative changes in head and 
plant efficiency in the long-range load allocation problem. 

Work is progressing on both the long and short range problems 
at the Bonneville Power Administration. A series of unpublished 
memoranda give the results of some of their work to date. Hoard 
(26) has given a clear description of factors entering the long 
range problem, together with an analysis of various modes of 
operating Grand Coulee, and an outline of the development program 
undertaken at B.P.A. (Whitbeck (6£) presented a report on the 
mathematical classification of the problem, identifying it with 
the calculus of variations. A mathematical formulation and solu¬ 
tion of the short range problem have been evolved by Hoard, Whit¬ 
beck and Wiener (Weiner and Whitbeck, 66, 67), which allocates 
load so that the sum of generations equals system load plus an 
estimate of system losses, with each hydro plant 



- 17 - 

constrained to deliver a specified average output for a short 
interval. As yet, this method causes violation of operating 
limitations because they have not been included in the mathe¬ 
matical formulation. There is, however, a strong likelihood 
that the use of fictitious costs as functions of limitation- 
violations will remedy this. A contribution to this develop¬ 
ment of a solution for the short range problem is presented 
in Appendix D, which does include limitations on plant 
capacities. 



Chapter III 

THE LONG-RANGE WATER UTILIZATION PROBLEM 

3.1 Nature of the Long Range Problem 

The determination of operations for one week at a time 
is not adequate when cumulative effects of weekly operations 
can materially influence plant efficiencies. This is the case 
in systems with large amounts of storage drawdown, where var¬ 
iations in head over relatively long periods of time result in 
changes in the amount of energy yield from a given flow through 
the turbines. This appreciably changes the character of the 
problem. This dependence of future plant efficiencies on 
prior operation requires that an optimization be sought over 
a sufficiently long future time interval. What is needed is 
not a set of load allocations for a given system load at any 
one time, but rather sets of time functions of generations 
over a long future time interval corresponding to the time 
functions of total system load. The dependence of plant 
efficiencies on elevation changes means that the principal 
variables are no longer plant generations but rather the specific 
combinations of elevation and discharge that determine these 
generations. These in turn require explicit use of natural 
stream flow and load predictions for this future time inter¬ 
val. The net head and discharge are further dependent on the 
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characteristics of the storage basin at each plant, and also 
on the characteristics of the tail water rise at each plant. 

In addition, with several widely separated plants on the same 
stream, the time of flow for water to flow between these 
plants must be taken into account. Thus a host of more ele¬ 
mentary system characteristics and functions must be introduced 
into this problem. 

Again, the shape of the cost curve for thermal or other 
non-hydro energy sources is of importance. The cost of such 
energy, in general, rises faster than linearily with increasing 
generation. In some cases the duration of such energy demand 
may also influence the effective costs of the hydro deficiency. 
Such would be the case, for example, if loads have to be inter¬ 
rupted because of inadequate thermal replacement for hydro 
deficiency. In each application of a long-range optimization, 
therefore, careful consideration must be given to the determin¬ 
ation of an effective cost for the replacement of hydro defi¬ 
ciency. The optimization, then, would seek to reduce the in¬ 
tegrated cost of operation over the time interval considered. 

3.2 Seasonal Variations 

The long-range problem is characterized by seasonal shifts 
of load requirements and non-coinciding seasonal variations in 
stream flows and hydro plant capabilities. Usually the peak 
in the over-all system load requirements and the lowest values 
of total natural flow energy occur in mid-winter. This com¬ 
bination often results in the necessity for appreciable storage 
drawdown and thermal use during this period. Such storage use 
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naturally reduces the energy yield from subsequent natural 
stream flow because of the reduction in plant head. This 
condition prevails until the reservoir is refilled, so that 
drawdown early in the season may cause a cumulative loss of 
energy of considerable magnitude. On the other hand, fail¬ 
ure to have a storage reservoir sufficiently drawn down at 
the beginning of the high-flow season results in spillage of 
water that otherwise could have been utilized for generation. 
Normally more than one storage reservoir must be drawn down 
simultaneously in order to meet the load demand or to"create 
a hole" for the spring flood. Proper coordination of draw¬ 
down is necessary in order to obtain best system operation, 
taking into account hydraulic connections at various storage 

dams and a time delay of water flow between dams. 

3.3 Accuracy of Flow Forecasts 

Of fundamental importance in the determination of long- 
range optimization is the accuracy of forecasts of stream 
flows and loads. Any one optimization can be only as accurate 
as the stream forecasts. It is imperative that such optimiza¬ 
tions be kept continuously or periodically up to date as often 
as fresh information on stream flow predictions becomes avail¬ 
able. 

Flow records for a long time in the past and ground-water- 
level measurements enable flow forecasters to make predictions 
that flows will exceed a certain minimum, with practically a 
100% probability. Various techniques can be used to predict 
the probabilities of the flow exceeding levels above this mini- 
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mum (Light and Kohler, 37; Stanley and Kennedy, 53; Bernard, 7). 
Whether operation is to be based on the expectation of just ex¬ 
ceeding the minimum flow, or on a lesser probability of exceeding 
some higher flow, is today a management decision. The use of 
conservative estimates, and consequent loss of energy through 
spillage, must be charged to insurance. Perhaps the greatest 
deficiency of the work presented here is that the various prob¬ 
abilities of flows are not used explicitly in the optimization. 
This eventually should be done, particularly when further de¬ 
velopments in the science of stream flow forecasting make avail¬ 
able stream flow probability distributions, as functions of 

time. However, in the following, it is assumed that only a 
"best” estimate of each flow is available, or that solutions 
may be carried out for more than one set of flows, such as 
median, minimum, and probable flows. 

3.4 Project and Operating Limitations 

Of great effect on the mathematical formulation of the 
long-range optimization are the project and operating limita¬ 
tions. The principal project limitations are the maximum dis¬ 
charge that can pass through the turbines at maximum gate open¬ 
ing as a function of the net head at each plant, and limitations 
on the minimum drawdown elevation because of the location of 
intakes or the limits of the storage basin. A less common proj¬ 
ect limitation involves channel restrictions which limit the 
flow of water as a function of water elevation. Operating 

%ork has been started on this feature by Mr. J. Little in a 
companion Ph.D. thesis. 
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limitations include minimum plant discharges for navigation 
or fish-life purposes, and maximum storage elevations as 
dictated by flood prospects. Many of these limitations are 
inequalities in form and present difficulties in the setting 
up of the mathematical models. 

3.5 Multiple Benefits 

A further generalization of the mathematical formulation 
may be necessary when the optimization is to be performed with 
respect to more than one system benefit. This is the case, 
for example, when irrigation water is required at times of no 
spillage. Then best operation would require economic and 
social balance between the benefits of power generation and 
irrigation. However, this is simply handled if it is de¬ 
cided to demand a specified amount of water for irrigation as 
a function of time and to then perform the optimization with 
respect to power generation having this demand as a constraint. 
Thus, in the following, it is assumed that irrigation has al¬ 
ready been accounted for in the prediction of natural stream 
flow availability. 

3.6 Length of Future Time Interval Considered 

A decision must also be made on the length of the future 
time interval to be used in the long-range optimization. To 
a considerable extent the economic operation will depend on 
requirements set down for the potential hydro energy remaining 
in the storage basins at the end of this time interval. This 
in turn must depend on a management decided risk to be taken 
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on predicted water resources beyond the end of this time 
interval. The selection of the end or the beginning of the 
spillage season as a future boundary has the advantage of 
simply requiring that all dams that spill be full at the end 
of the time interval. However, in those dams which do not 

spill, it is necessary to decide whether a net increase or 
decrease in storage is desirable during the future time in¬ 
terval. In a fully developed system where a large number of 
plants do not spill during any portion of the year, time in¬ 
tervals greater than one year may be advantageous. 

There is some question concerning the effect of increas¬ 
ing uncertainties in flow prediction, as the length of the 
future time interval is increased. A point of diminishing 
returns may be reached in some cases as the length is in¬ 
creased. However, the nature of the problem seems to require 
that the interval include at least the forthcoming refill 
season, because the ability to refill using water that other¬ 
wise would be spilled is the justification for draw-down. 



Chapter IV 

THE ROLE OF THE COMPUTER 
IK THE OPTIMIZATION OF SYSTEM PERFORMANCE 

4.1 Introduction 

In electric power systems, several types of supervisory- 
control are exercised to achieve economic operation. It is 
the purpose of this chapter to consider the nature of such 
supervision, from a computer-control viewpoint. Although the 
ingredients of control systems vary in different electric 
systems, depending chiefly on the nature of the energy re¬ 
sources used, certain basic concepts apply to each, and are 
discussed. 

4.2 The Process Function 

The process we are concerned with is the conversion of 
energy into electrical form. The cost of operation is depend¬ 
ent upon a great many factors, but only a few of these ordinar¬ 
ily are variable and need be considered in the optimization of 
performance. The functional dependence between the performance 
of the system and these pertinent characteristics of the actu¬ 
ating inputs alone determines the "process function" as far as 
the optimizing controller is concerned. This relationship is 
sketched in figure 4.1. Parameters of this functional relation¬ 
ship are the load allocations, which can be adjusted by a con¬ 
troller so as to achieve an optimization of performance. The 
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process function and its actuating inputs are all a computer- 
controller need know about the physical system being controlled. 

4.3 Thermal System Controller 

Simplified block diagrams of system economic controls are 
shown in figures 4.2, 4.3, and 4.4. In figure 4.2 is the block 
diagram of a purely thermal system control. The performance 
variables in this case are fuel costs and the operation of 
plants within their capabilities. The pertinent input char¬ 
acteristic is fuel prices and the adjustable parameters are 
load allocations. In the case of urban thermal systems, the 
process function must contain information on fuel cost char¬ 
acteristics of each thermal plant. In the case of interconnected 
thermal systems, the process function must, in addition, contain 
information on transmission losses as functions of load alloca¬ 
tions. The controller provides instantaneous load allocations 
in accordance with the instantaneous total system load. 

4*4 Short Range Hydro-Thermal System Controller 

In Fig, 4.3 is shown a block diagram for the economic con¬ 
trol of a hydro-thermal system having run of the river plants 
and pondage plants with no appreciable cumulative variation in 
plant head and efficiencies. The controller dictates load allo¬ 
cation based on an optimization over a short future time interval. 
Instantaneous load allocations are given in accordance with the 
instantaneous total system load, and knowledge of the weekly 
load pattern. The optimization by the short-range controller 
is in turn dependent on the specification of the weekly allot- 
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ment of water for each hydro plant. This is determined at a 
higher level evaluation of past flow records, extrapolated 
into the future time interval. 

4.5 Long Range Hydro-Thermal System Controller 

In Fig. 4.4 is shown the block diagram of the economic 
control of the system where some plants experience cumulative 
changes in heads and plant efficiencies during a long time 
interval. Storage elevations are added as pertinent inputs 
to the process function. Benefits other than purely electri¬ 
cal generation must often also be considered as performance 
variables. The short-range controller in this case is in¬ 
fluenced by the weekly water allotment dictated by a higher- 
level long-range controller. The long-range controller in 
turn receives as inputs predictions of flows and loads over 
the long-range future time interval. This process function 
includes, in addition to the cost characteristics of the 
thermal plants and transmission loss relations, information 
on storage basin characteristics, variation of hydro-plant 
efficiency with net head and discharge, and the tail water 
characteristics of the hydro-plants. 

The long-range hydro-thermal system-controller is the 
special concern of this treatise. Therefore it will be fur¬ 
ther elaborated on. A pictorial sketch of such a control 
system, in a simple case involving two hydro plants and one 
source of thermal power, is shown in figure 4.5. In this 
integrated system, there are a number of semi-independent 
control systems superimposed upon one another so that each 
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in turn forms a part of a higher level control system. Within 
each plant there are many regulators to maintain a given oper¬ 
ation. Above these are the supervisory control for allocation 
of load among units in the plant. Each plant then becomes 
part of a control system under the direction of the dispatching 
office which allocates hourly loads on a short range basis. 

This unit in turn is part of a higher level of long range 
supervisory control which evaluates long-term operation and 

determines the amount of water to be made available to the 
dispatcher for the week’s operation. 

At each level, the supervision is concerned not with the 
details of lower level controls but only with overall character¬ 
istics. Moreover, the time scale of operations changes with the 
level of control so that at each level operation is based on 
time-averaged behavior of lower-level controls. Thus, at the 
dispatching level, an hourly average load is assigned; and at 
the next level the weekly allocation of water is made. The 
latter is dependent on the proper functioning of the former, 
because the average efficiencies of each plant depend on the 
dispatching practice used. Hence, the long-range water allo¬ 
cation is dependent on the weekly performance of the dispatch¬ 
er's short range allocation, and in particular on the weekly 
average energy yield per volume of water consumed at each hydro 
plant and the weekly average cost per mw hour at each thermal 
plant. 

At a given time, the state of the system water resources 
is in part a consequence of past orders from the controller. 
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In determining the weekly water allocations over a future time 
interval, the controller must take into account the present 
state of the system, the average operating characteristics of 
hydro plants and thermal plants as controlled by the dispatch¬ 
ing center, the distribution system, and predictions of load 
and stream flow conditions for the future time interval. 

4.6 Effective Cost of Operation 

The objectives of each controller must be clearly defined 
in terms of the minimization of an effective cost of operation. 
To obtain a single index of cost of operation, the concept of 
cost must be generalized to take into account all performance 
variables. The effective cost must be an aggregate of all 
factors considered to be of economic or social significance 
weighted according to their relative importance. To obtain a 
quantitative measure of cost of operation, an equivalent dollar 
cost must be assigned to each of these factors. In particular, 
the undesirability of the violation of project and operating 
limitations can be taken into account by the assignment of ef¬ 
fective penalty costs when such violations occur. These costs 
as functions of the degree of limitation violation must be 
designed so that only operations with an arbitrarily small 
amount of violation of limitation will be economical. The 
effective cost per hour then is a summation of terms, one of 
which is the dollars per hour fuel cost of the thermal plants. 

4.7 The Performance Surface 

The cost index is functionally related to all of the 
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process-function inputs, that is, the pertinent characteristics 
of the actuating inputs, and the adjustable parameters. The 
objective of the controller is to minimize the cost index by 
proper adjustments of the adjustable parameters. A geometrical 
interpretation of this minimization enables one to visualize 
the problem. In three-dimensional space, for example, a sur¬ 
face is defined by a function, z(x,y), such that the values 
of x and y as points along two axes locate a point on the sur¬ 
face, whose projection onto the third axis is equal to z. A 
fixed value of y determines a curved line on this surface. 

The curve of z versus x with fixed y will have certain minimum 
values of z, and x can be adjusted to seek these minima. This 
same concept can be extended to more dimensions. One coordinate 
is, of course, the cost index which is to be minimized. The 

other coordinates are the variable relevant characteristics of 
the actuating inputs and the adjustable parameters. For a 
given set of actuating inputs, the parameters are to be adjusted 
so as to minimize the cost index. As the actuating inputs 
change, operation takes place in a different region of the sur¬ 
face, so that a revised set of parameter adjustments is needed 
to minimize the cost index. In many physical problems the 
shape of the minimum is relatively broad and only one minimum 
exists within the limitations of the problem. 

A slightly different interpretation is to consider only 
the adjustable parameters as coordinates in this multi-dimen¬ 
sional space so that the relevant characteristics of the actu¬ 
ating inputs are included in the functional character of the 
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surface itself. The surface then changes as the relevant 
characteristics of the actuating inputs change, but for a 
given set of actuating inputs, the surface is stationary. 

The parameters are then to be adjusted to seek the actual 
minimum in this surface with respect to the cost index 
ordinate, A single point on either surface determines a 

mode of operation, 

4.# Computer Procedures 

Consider a mode of operation that differs somewhat from 
the optimum. In the surface having only adjustable parameters 
as coordinates, the operating point is somewhat removed from 
the surface minimum. Two courses of action are then open to 
the controller. It can determine the location of the minimum 
for the given set of actuating inputs. Then, the parameters 
can be adjusted in one step so as to move the operating point 
to the minimum. Alternatively, the controller can determine 
the path to be followed on this surface toward the minimum 
and can adjust the parameters in steps so as to move step-wise 
towards the minimum. Computers following these two procedures 
will be referred to as destination-finding computers and path¬ 
finding computers respectively. 

In the case of load allocation among thermal plants, the 
necessary conditions for the gradient of the surface to be 
zero can be solved*, in order to obtain the load allocations 
for different system loads. The results of these computations 
can be stored in a variety of forms. The function of the 
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controller may then be simply to obtain from storage the load 
allocations corresponding to the current system load. This 
is clearly a "destination-finding" type of controller. This 
obtaining the solution from storage has in fact been mechanized 
in the case of allocating load among units within a plant 
(Purcell, 44; Purcell and Powell, 45). Conceivably, it can 
also be mechanized on the system level. 

In the case of purely thermal systems, with fixed fuel 
prices, the only input that must be carried to storage is 
total system load. In systems containing hydro plants, how¬ 
ever, the short range controller inputs also include water 
allotments for each hydro plant for a short future time inter¬ 
val. While it is still conceivable to store the load alloca¬ 
tions for all expected combinations of the controller inputs, 
the multivariable storage requirements would be very great 
with more than a few hydro plants. It appears that under cer¬ 
tain assumptions, in this case, it is possible to obtain an 
analytic solution (to the equations setting the gradient equal 
to zero) in terms of the controller inputs (Weiner and Whitbeck, 
66 and 67). Hence, a simpler controller would store the solu¬ 
tion in terms of this functional relationship. The controller 
would then operate by computing this function for the current 
controller inputs. This still would clearly be a destination 
finding controller. 

In the long range optimization of hydro-thermal system 
operation, it appears that analytic expressions for the solu¬ 
tion to the equations setting the gradient equal to zero, can 
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be obtained only with unrealistic assumptions and unjustifiable 
simplifications. This is true primarily because of the addi¬ 
tional unilateral project and operating limitations in hydro 
plants, that must be included to obtain realistic operation. 
Hence neither the storage of complete solutions or the storage 
of analytic expressions for the solutions appear feasible in 
this case. What can be stored are the determining equations 
themselves; that is, the expressions for the gradient of the 
surface. The long range controller must, therefore, store 
these determining equations, in either analog or numerical 

form, and re-execute their solution for each new set of con¬ 
troller inputs. 

Whether the gradient is evaluated to determine the path 
to be taken, or whether the minimum is obtained by setting the 
gradient equal to zero, the possibility of more than one mini¬ 
mum must be contended with. The destination-finding computer 

using the path of steepest descent is best used to periodically 
correct for small deviations from the best operation. This is 
the case, for example, when periodic corrections are made dur¬ 
ing a season as revised predictions of stream flows become 
available. In general, with this method, it is necessary that 
a first approximation to the best operation be available which 
will set operation within the ’’bowl” of the correct minimum. 
Similarly, if the equations setting the gradient equal to zero 
are solved, more than one solution may result. The one with 
least cost and yet physically realizable must then be selected. 
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4.9 Path of Steepest Descent 

When moving the operating point toward the minimum with 
a path-finding computer, one intuitively would make changes 
in each adjustable parameter in proportion to the reduction 
in the cost index achieved by a small change in each adjust¬ 
able parameter. That this would move the operating point 
along the path of steepest descent is demonstrated in Appen¬ 
dix A. The magnitudes of the components of the vector: 



where * the adjustable parameter 
M ± z unit vectors 

give the relative improvements in performance to be gained 
by a change in each of the adjustable parameters. The alge¬ 
braic sign of each component determines the direction in 
which it would be advantageous to adjust each parameter. 
(Appendix A). 



Chapter V 

DEVELOPMENT OF THE DETERMINING EQUATIONS 
FOR THE LONG RANGE SOLUTION OF THE THREE PLANT SYSTEM 

5.1 Introduction 

The basic ingredients of a solution to the long range 
problem are (a) the mathematical formulation of the problem 
and (b) the mathematical development to the point where the 
solution rests only on the evaluation or solution of a 
specified set of equations. These equations, which determine 
economic operation, are referred to as the ’’determining equa¬ 
tions”. 

To avoid unnecessary complexity, while retaining the 
principle elements in the problem, the determining equations 
are first found for the simple three plant system shown in 
figure 4.5. They are then generalized to apply to larger 
systems. 

5.2 The Cost Index 

The effective hourly cost of operation is defined as the 
hourly cost of replacing the hydro-deficiency plus penalty 
costs, depending on storage elevations and plant discharges, 
for violations of project or operating limitations. Thus: 

t * 

$/hr * $/hr (g d ) ♦ £ P iy + £ P id (5-1) 
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Sketches of penalty functions, which must be designed to make 
even a small violation of limitations uneconomical, are shown 
in figure 5.1. (The maximum generation or discharge is, of 
course, a function of the net head.) The objective of econ¬ 
omy control is then the minimization of the integrated cost 
over a future time interval, that is to minimize 


$ = 



r (t)dt . 


(5-2) 


5.3 Plant Characteristics 

The cost index is dependent on various hydro and thermal 
plant characteristics. These, sketched in figures 5.2 through 
5.11, are: 

(a) tailwater elevation vs total plant discharge 

(b) plant efficiency vs turbine discharge with net head 
as a parameter 

(c) spillage vs storage elevation above crest 

(d) storage vs storage elevation 

(e) maximum plant generation (or discharge) vs net head 

(f) incremental cost for replacement of hydro deficiency 
vs hydro deficiency. 

These are largely based on actual data, supplied by B.P.A., 
for specific installations. In those cases where inadequate 
data was on hand, attempts were made to make the functions 
flexible enough to characterize a variety of plants. The plant 
efficiencies, in particular, were so formed. 
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It was explained in section 4.5 that since the long range 
solution seeks only weekly water allocations to each hydro 
plant (and the corresponding energy required from thermal 
sources) it is concerned only with average plant character¬ 
istics for that period. Instantaneous hydro plant efficiencies 
are shown in figures 2.2 and 2.3 as functions of plant output. 
However, for the long range optimization, an averaged set of 
curves consistent with operating practice should be used. For 
a given number of machines on the line in a hydro plant, at a 
given head, figures 2.2 and 2.3 show there is one discharge 
which will give highest plant efficiency. Good operating 
practice is to operate at or near these points. The flatness 
of these peaks is dependent on the type and number of machines 
used. Adjustable blade turbines, for example, have relatively 
flat characteristics. Hoard has pointed out (26) that with 
many units on the line, operation can always be had at very 
nearly peak efficiency. With 1065 raw load at Grand Coulee, 
he illustrates, operation of 10 units above their points of 
maximum efficiency or operation of 11 units below their points 
of maximum efficiency requires the loading of each unit only 
5 per cent from their most efficient load, with a reduction in 
plant efficiency of only 0.12 per cent. In any case, it is 
reasonable to assume that the short range load-dispatching pro¬ 
cedure will usually load all plants near points of peak effi¬ 
ciency. Loadings at different hours of the day may be at dif¬ 
ferent peaks, so that an average efficiency between these 
points may result. Therefore, the average plant efficiency at 
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a given head is characterized as a function of discharge by 

a smooth curve, passing through the regions of highest effi¬ 

ciencies in figures 2.2 and 2.3. In order to allow consider¬ 
able flexibility, which might be needed in plants of widely 
different characters, these curves are expressed as: 

. »li [hi(t) . h io _ m ig D i(t ) . n ig D l(t ) 2 ] ( 5 - 3 ) 

Plots of these curves, for cases with somewhat exagerated 
variation of efficiency with discharge, are shown in figures 
5.4 and 5.5. 

The replacement of hydro-deficiency is provided, in the 

3 plant problem, from a thermal plant. The cost of such energy 

has three major components: one which is independent of the 
load (such as overhead and capitalization), one which is pro¬ 
portional to load, and one which varies as the square of the 
load (due in part to losses, reduction of efficiencies at higher 
loadings, and limited capacities of equipment). The cost char¬ 
acteristic used for this problem, figure 5.12, corresponds to 
that of a relatively large, efficient, thermal station, taken 
from the plant data given by George Page and Ward (22). 

5.4 Establishing the Functional Character of the Cost Index 

In what follows, storage and rate of change of storage 
are treated as the fundamental variables. Alternatively, stor¬ 
age elevation and its rate of change could be used. However, 
in different plants a foot change in elevation can involve 
vastly different changes in amounts of water, depending on the 
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storage basin characteristics. Hence storage is considered 
more basic. At any rate, the two are related by storage 
characteristics, figures 5.3 and 5.9, and by the incremental 
storage in: 

s (y(t)J jf (t) * Stt) (5-4) 

The three plant system of figure 4.5 is now used as a pilot 
system to establish basic relations. The storage at dam F, at 
a future time t is: 

S,£t) - + /[ROO -D.Cx) <kx 

° (5-5) 

Continuity requires that at each plant the outflow equal the 
inflow minus the rate of change of storage, i.e.: 

D,ct) + t) - RCt)- 5,Ct) (5-6) 

The inflow to plant 2 (see fig. 4.5) includes the output from 
plant 1 at a previous time plus the tributary inflow (or run¬ 
off between plants). Hence, as above: 

S 2 (t) = $ 2 Ct=o) + J [ ff*(x) + D» (x-Vz) + 

- D a oo-JU *)]dx (5 ” 7) 

and also: 

+ 6ft)+D ,- S'W (J.8) 

Then, substituting (5-6) in (5-3): 

0 z ct> * FTmCt) + £ Ct-Um) - W ( 5 _ 9 ) 

Consider g-^t). From figure 5.4, it depends on h 1 and D^. 
By definition, net head is the difference between storage eleva¬ 
tion and tailwater elevation. Hence g^(t) depends on y ^(t), 
y 1T (t) and D]_(t). From figure 5.2, y, T is dependent on 
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Q]_( t), equal to turbine discharge plus spillage. The latter 
is, from figure 5.6, dependent on storage elevation alone. 

Hence we know that g^(t) depends only on y,(t) and D^(t). 

From equation 5-6, turbine discharge is dependent on natural 
inflow, rate of change of storage, and spillage. Again, the 
latter is dependent only on Jf, (t). Hence g^tt) depends only 
on F^Jt) and S, (t). These are fundamental quantities. 

Since if, (t) is directly related to S,(t) by figure 5.&, we 
establish the functional character of g-^(t) as: 

Sl(t) * g 1 ( F 1 (t), S 1 ( t), SjJt)) (5-10) 

In like manner, it is easily shown that the functional char¬ 
acter of g 2 (t) is: 

g£(t) s §2 ^ ^i (b~ T 12 ) i F-^t i) i S^(t- T^)* ^2 (t)» 6 2 (^ ^ 

(5-11) 

By definition, the hydro deficiency is the difference between 
system load requirements and the total hydro generation. In 
this case, therefore: 

g d (t) * L(t) - g 1 (t) - g 2 U) (5-12) 

Consequently, using both (5-10) and (5-11), the hydro deficiency 
is seen to depend on: 

g d (t) = g d (ut), F x (t), F 1 (t- r 12 ), F 12 (t),s 1 (t)s 2 (t) 
S x (t), S x (t- T 12 ).S 2 (t)) (5-13) 

Since the penalty costs, sketched in figure 5.1, are likewise 
functions of turbine discharge and storage elevations, it fol¬ 
lows from (5-1) that the effective $/hr cost has the same 
functional character as in (5-13). 



- 46 - 


5.5 Setting Ug Changes in the Cost Index by the Calculus 
of Variations 

The objective is to determine the curves S 1 (t) and S 2 (t) 
for the time interval chosen which will minimize the cost in¬ 
dex defined in 5.2 as: 


$ = 



r(t)dt 


(5-2) 


No generality is lost in specifying that S 1 (t) and S 2 (t) are 
continuous and have continuous first and second derivatives* 
Then, since all system characteristics are at least sectionally 
continuous, #/hr(t) will be sectionally continuous. Hence, the 
conditions of admission for problems in the calculus of varia¬ 
tions are fulfilled (Courant, 12). 

In practice, the optimization is to be carried out over a 
future time interval. Hence, for t^-0 , the storages are all 
fixed, and no variations of the storages are allowed. Also, 
the end of the time interval must be interpreted in terms of 
the physics of the problem. If the operation of plant 1 is 
specified up to time T- f^ 2 , the effects of this operation 
will be felt at plant 2 up till time T. Plants further down¬ 
stream would be affected at still later times. To consider 
the delayed consequences of operation, then, the time interval 
must be extended for the downstream plants. The upper limit 
on 5-2 is to be so interpreted. Hence, in modifying a given 
set of storage curves in order to improve operation, we are 


constrained so that: 


SS.(t<o) = iS,(tio) = o 


( 5 - 14 ) 



- 47 - 

£ St (t-^T) =* 0 ( 5 - 15 ) 

<rs,(oT-r,i) = o (5 _ 16) 

The functional character of the integrand of (5-2) has been 
established in 5.4 as: 

$/hr (t) = s,(t), 

S t (t)t SiCt'fci), Sj(t)J (5-17) 

In the usual manner of setting up such problems (Courant, 12), 
let the actual storage curves be set equal to the optimum 
curves plus an "error function”, i.e., let: 

S, (t) = s..(t) + * <p(t) (5 _ 1S) 

& (*) - S,.(t)+ tr w C*> (5 " 19) 

S 1Q (t) and S 2Q (t) are the curves sought, which give the mini¬ 
mum value of the cost index, (p (t) and ((J (t) are fixed but 
undetermined curves which are constrained by equations (5-14, 
15, 16) so that 

t>(t^o) * 4) (t^O) = o (5-20) 

0 (t ^T~Vz) — 0 ( 5 - 21 ) 

V(t»T ) 9 0 (5-22) 

zt and v are variables whose magnitudes determine the differ¬ 
ences between the actual and the sought-for curves. It follows 
that changes of M and v can bring each storage curve closer 
to the optimum curves. 

In (5-17), the only terms dependent on m are S 1 (t), 

S x (t )f and S^(t— T ^2)• Hence: 

trlSJ dmr. 1 , JtlttL ggS' i 9S/kr- w 

P$,ct-W 

sMk 


9 m 


~h 


9$,Ct) 3 m ' r 3S,(t) 3 m 

- as ♦»+as ♦<*>+ & 


(5-23) 
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Now: 


it ( 5 - 24 ) 

Substituting (5-23) in (5-24), the second and third terms are 
integrated by parts. For the second term: 


(«*$&*■ 


(5-25) 


Similarly, the third term of (5-24) gives, on integration: 

At 




(5-26) 


Using the constraints of (5-20) and (5-21), it is seen that 
the integrated parts of (5-25) and (5-26) are zero. The in¬ 
tegration of ( 5 - 23 ) thus becomes: 

- f T tb(t) f - - 4 - 2&221 it 

J* ss.ct) J 

it 


rV rt n [ j_MlkCiL 
- J f(t w [ Jt PStt- 77 ') 


(5-27) 


In like manner, since in (5-17) only S 2 (t) and S 2 (t) depend 
on v : 


lilkM = ymlMdil + juffl 4^r 


(5-23) 


The second term of (5-23) is likewise integrated by parts. 

This using (5-20) and (5-22) to eliminate the integrated part, 
one gets: 


dt =-{Vt 


Jt 


(5-29) 
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Again, 


f Jlkm Jt - JJL 

4 Jv ~ 2v 


Therefore, using (5-29) in (5-28) in (5-30) 


(5-30) 


(5-3D 


The total change in the cost index is: 

+ 0 Sr 

Substituting (5-27) and (5-31) in (5-32) one gets: 

= Sm ^ ^iJhrC H „ A- S$/hr dt 


B<4, i vw l 7 &w 

- Sm j [jt pfft-m)] 


(5-32) 


(5-33) 




Because and are independent of t, they can be brought 

inside the integrals. Furthermore: 


4>(9 = -jjf* 1 Sm = / & (TO 

Sm = - SJ ^ n ^ Sm = SS,(t-K) (5 _ 34) 

<rv = - <r^, cd 

d i 

Equation (5-33) then gives: n 

. r T e c g) \iHkiS \_ jt 

In order to combine the first two integrals in (5-35), with 
£ S-^(t) as a common factor, the integrand of the second inte¬ 
gral must be shifted in time an amount Tdays. The limits 


(5-35) 
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of integration must be changed accordingly. Thus: 

- -r 

•-fjj 

As the problem has been set up, however, time t <0 is past 
time for which S S-jJt) is constrained to be zero. Hence the 
lower limit of integration in (5-36) may be set at zero. Also, 
using (5-15)i the upper limit of integration in the first in¬ 
tegral of (5-35) may be set at T- T^. Equation (5-35) is 
then rewritten as: *>// . 

r T ' r " r s&krd) J Jill*® - A- 

ft = l - It JS.O) at 3 s.ct) J 




t it (5-37) 


This equation relates the storage corrections, made throughout 
the interval considered, to their net effect on the system 
cost index. 


5.6 Role of the Gradient 

As described in 4.#, the determining equations may be 
used in either of two ways. They may be used to express the 
necessary and sufficient conditions for a solution to exist. 
That is, in terms of the surface of 4.7, they may express the 
conditions for a stationary point on this surface. These 
conditions are that all components of the surface gradient be 
zero. Alternatively, the determining equation may be used to 


c/t 


*A stationary point is a maximum, minimum, or saddle point 
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determine a path to be followed towards a solution, or towards 
a minimum in the surface of 4.7. When the path of steepest 
descent is to be followed, the determining equations should 
give the components of the gradient of the surface at the 
operating point. Thus, for either procedure, we need expres¬ 
sions for the gradient of the surface. 

We seek to optimize operation over a future time inter¬ 
val. Hence we seek time functions of storage throughout this 
interval. These continuous curves contain a finite amount of 
information, depending on the highest frequency components in 
these curves. Hence, for these relatively smooth curves, 
knowledge of the storages at a finite number of points (to 
within a reasonable tolerance) adequately characterizes them. 
These ordinates of the storage curves are, then, the variable 
parameters in the optimization. The components of the surface 
gradient are then nothing more than the changes in system cost 
index per unit change of each of these ordinates. 

5.7 Proof of the Determining Equations 

From equation (5-37) the components of the gradient and 
hence the determining equations can be derived. This was in¬ 
dicated by Courant (13) who observed the correspondence, in 
calculus of variations problems, between the gradient of a sur¬ 
face and the Euler expression, typified by the bracket in the 
second integral of (5-37). Two derivations of the determining 
equations from (5-37) are given in this treatise. Because of 
their length and the reasonableness of the results, these 
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proofs are relegated to the appendices. The first, presented 
in Appendix B, involves a geometrical interpretation in multi¬ 
dimensional space using the analogy to a surface gradient. 

The second, presented in Appendix C, does not involve the 
assumption of a finite number of storage ordinates, or the 
gradient concept, and relies on bounding the changes in cost 
index using Schwartz’s inequality. 

Both proofs involve the identification of the gradient 
with these changes in the storage curves (throughout the time 
interval considered) which produce the maximum reduction in 
the cost index per unit of ’’integrated storage change squared”. 
In the more general proof of Appendix C, the latter is defined 




(5-3$) 


In the geometrical interpretation involving a finite number of 
ordinates, presented in Appendix B, the corresponding ’’summed 


storage-change-squared" is used: 

URl 1 -if 


M. ordinates/plant 


,«/ 


n plants 


(5-39) 


Here | S R| is the size of step taken in a performance surface-s¬ 
having the storage ordinates as orthogonal coordinates. Thus 
the maximizing of the ratio: reduction of is synonymous 

with taking the path of steepest descent and proceeding opposite 


1 


section 


4.7 
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to the direction of the surface gradient. The final results 


are as follows: 

a) In order to obtain the largest reduction in the cost index 
per unit of integrated storage-change-squared, i.e., to 
move along the path of steepest descent, storage changes 
throughout the time interval considered should be made 
according to the following formulas: 


/ S. (t) - 


< 


3$kr(t) _ J$/hrCt) j 3tlKi*$ 

dt Jtxf) I' 5 "'* 0 ’ 


S $,(t)= - 



3$/kr(i) 

3 Si C t) 


dt JSzCV 


(5-41) 


where »( is a positive constant determining the size of step 
taken. 

b) The necessary and sufficient conditions for a solution to 
be stationary is that the above equations be equal to zero 
throughout the time interval considered. 

These expressions are further discussed and developed in the 
following chapter. 



Chapter VI 

GENERALIZATION AND FURTHER EVOLUTION 
OF DETERMINING EQUATIONS 

6.1 The Determining Equations for Larger Systems 

The results of chapter V can be generalized to apply to 
an arbitrarily large system. The basic objective is kept 
the same, that of reducing the cost of hydro-deficiency re¬ 
placement. It is presumed that a given hydro deficiency 
could be divided up among the non-hydro sources so as to 
minimize its cost. If thermal plants supply the deficiency, 
the incremental loading procedures outlined in chapter II 
could be used. There exists, then, in general, a curve of 
minimum cost of hydro replacement as a function of mw hydro 
deficiency, as in chapter V. 

There remains to be considered the effects of additional 
hydro plants on the same and different streams. These addi¬ 
tions add nothing fundamentally new to the problem, and their 
effects can be ascertained from the interpretation of com¬ 
ponents in equations (5-40) and (5-41). Naturally, changing 
the discharge of an upstream plant affects operation of the 
downstream plants. This is reflected in the equation for 
storage change, (5-40), whose major bracket has, from eq. 
(5-37), the dimensions of $>/hr/ksfd. The influence of / S,(t) 
on operating costs at time (t), due to affected output of 
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plant 1, is given by the first two terms of the bracket in 
(5-40): 

/ f/kr(t) _ 4 §)hr (1) _ 4 

jS.Lt) ~ JS.U) dt 3S,Ct) (6-i) 

The influence of S S,(t) on operating costs at a later time, 
(t*r l2 ), due to affected output at the downstream plant is 
given by the third term of this bracket: 

/ t/hrCt +Tn) _ d ji/irCt+Vi) , 

S S. Ct) dt 3 s,Ct) l6 ' 2) 

Additional downstream plants would logically cause additional 
changes in operating costs at still later times, so that addi¬ 
tional terms of the type in (6-2) would be added to equation 
(5-40). In the three plant problems, the evaluation of (6-2) 
depends on the presence of the tributary between plants 1 and 
2, but this presence does not affect the form of (6-2). Sim¬ 
ilarly, the terms added to (5-40) to account for additional 
downstream plants have the form of (6-2) regardless of the 
presence of other joining streams. Regardless of whether the 
streams of a system are arranged in parallel or series- 
parallel, each plant finds itself in the position just 
described where changes in its output affect operation at a 
series of ’’downstream" plants. (Plants downstream from plant 
i are simply those through which water from plant i passes). 

A change in the operation of any plant must rest on an evalua¬ 
tion of the effects of such change on all downstream plants. 
Hence, the general expression for the correct change at each 
plant i is obtained by direct extension of (5-40) as: 
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SSi(t) = 

i = 1, 2, ....n hydro plants.* 

Here the summation must be interpreted to include the ith 
plant and all plants downstream from it. It is understood 
that the addition of other downstream plants involves de¬ 
layed consequences of upstream plant operation, so that the 
time interval under consideration must be extended for down¬ 
stream plants, as discussed in 5.5 . An expanded form of 
(6-3) is developed in the following. 


H/hr(ti J_ )$/kr(UHtj ) 

sSid) Z_ jt ap) J 

b*t J 


6.2 Analytic Expressions for the Partial Derivatives 

It is advisable, to reduce truncation and round off 
errors, to use analytic expressions instead of taking differ¬ 
ences to evaluate the partial derivatives in (6-3). These 
are now derived. 

a) "Basic Relations”. The following expressions are used to 
fit the characteristic curves presented in section 5.3. 
Tailwater Elevation: 

iJir = HiT(Q; = °) -+ ntirlDitJiJ + n ir [V;+J<] (6-4) 


Spillage: 


J { = Si + n iL S t 


Generation: 


gi = 7«' [ *jl ~ iliT ha ' Wig R “ n ('g V,' ] Di 

Storage Elevation: 

Hi - Hi. + S < + n <* S ‘ 2 4 °‘*» S ‘ 3 


(6-5) 


( 6 - 6 ) 


(6-7) 


* In the generalized equations, plant 1 is the furthest downstream. 
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The logical extension of the continuity equation (5-6) is 
given by: 

A (t) + J-(t) =£ - IT S, (t'Tqi) (6. 8) 

«=t uTJ 

Where the summation includes plant i and all upstream plants, 
defined as those plants whose discharge passes through 
plant i. It follows that for a plant downstream from plant i 


n n 

D b (t+r ib )= £ S„ft+U-r.b) 

u-b 


(6-9) 


These show the fundamental dependence of turbine discharge 
on the basic variables 5; (t) and Sj (t) to be: 


HIM = _ j 1 ffl - _ r (t) 

3 S, ct) j & c t) - ^ 
£Ml!*) = o b#i 

P $; Ct) 

P ^ Cf) 


( 6 - 10 ) 

( 6 - 11 ) 

( 6 - 12 ) 


b) "Effects of Storage Changes on Plant Generations." Combin¬ 
ing equations (6-4) and (6-6), one gets: 


Si 


<.(t) =fgj BKO = 7,{ j.f«-WirfD. +^ W] 

Lsh 1 , ,, 

- hi . - m,« P, ft) - % D; ft)'} P; (0 


( 6 - 13 ) 


The dependence of g^(t) on the basic variables Si (t) and 
S t - (t) can now be found. Let x represent either of these 
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variables. Then: 


W _ a^d. u), JaiJJMl + r .p, 

- B ; (r.«X 4 fti) 

A(lK*U«) = f£](0 - 7,{(m, 4 fm ir ]p,.(t) +2 n,gD,ft)‘ 

+ 2n ; rRC*)[l!ft)+Jft)]l 


where: 


(6-15) 


B,' Ct) =• i; D; Ci)^ m iT ■+" ^ n,- T [ D,- (fc) +Ct)j | 


( 6 - 16 ) 


s i 

tr-(t) is defined by equation (6-13) 
u i 

Since Sb(* +T *«>) is independent of S jL (t) 

3 HkCtiTj _ }Jt(t*ri t \ _ 0 

P 2; Ct) P X,- Ct) 


and S.(t): 

l . 

(6-17) 


Therefore, one gets from (6-14) and (6-17) the following 
dependence of future generation on present changes in 
S i (t) or S ± (t): 


sejt+u) 

3 Kid'S 


= A k (J, 4 (Wb) 


(6-18) 


The expressions presented in 6.2a are now made use of. 
Letting X i (t) = S ± (t) and substituting (6-10) in (6-14) 
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^ iCi l ,7 J Ml 

s, a) 6 w&c*) 




(6-19) 


where: 


c-ct) 5 A CO+ B.-CO 


= Ml d) - ij r (Mj e D, ft)+2 fligRft) 1 

i . 


( 6 - 20 ) 


Using (6-11) and (6-1&): 

f>hCt+71 k) __ Q 

P S,' Ci) 

Using (6-12) and (6-1#): 


bi^ i 


( 6 - 21 ) 


jLMiSI = _ A h (p t m t ),JMS) (6 . 22) 

* S; tt) 

Thus equations (6-19), (6-21), and (6-22) determine the effects 
of storage change on plant generations. 

6.3 Final Form of Determining Equations 

The hourly effective cost of the enlarged system is de¬ 
fined similar to (5-1) as: 


$A w = $/k r (s/j -t £ p !s (s .(4 + £ p, 

i s < i sl 


V 


(6-23) 


where the summations are over all hydro plants. Also, by def¬ 
inition of hydro deficiency: 


S’ft) = LCD -X>Ct) 


(6-24) 


t =£ b 


From ( 6 - 24 ) and the physics of the problem, it is known that 

= o 

„ S; CO 

and ( 6 - 24 a) 

^ gd (t+ftt) _ _ 

d ^ co 


(6-24a) 
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Neglecting the penalty cost terms in (6-23) for the moment, 
and using ( 6 - 24 a), equation (6-3) becomes: 


S S;(t)=+« 




3 & Ct) 




(6-25) 


where 


y(gj*>)= 


A $lht (&) 

A & 


( 6 - 26 ) 


sthe incremental replacement cost. 
Using (6-19) and (6-22) in (6-25), one gets: 

/(?/«)[ 7,4^ D,(t) - 

I d ( 6-27 

~h y(^d(t+T?4 AbCt+^b) \ 

bsl 1 

Next, the penalty cost terms of (6-23) must be included in 
(6-3). The additional terms are: 



dfjuJM 


I A S; Ci) 


t A Fud) 3b;(t) 

dDi Ct) 5SfCt) 



Using (6-10) and (6-12), this becomes: 

- * | p is ' ct) - fio (t) 4 'ft) + < 6 - 29 ) 

Adding (6-29) to (6-27), the complete expression is, finally: 


^S ; ct) = ^^(«,(a)[v ; jlfjRco-oau'ft)] - /j;c# 

+ f>.;c*)i'a)+i-^ hfccwaJVtoi)- p k ;ct^») j 


( 6 - 30 ) 


i = 1, 2, ...n hydro plants 
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The indicated summation includes the ith plant and all down¬ 
stream plants. One such equation exists for each hydro plant 
in the system. This completes the derivation of the general¬ 
ized "determining equations", which determine the most econom¬ 
ical operation of the system as described in section 5.6. 
Equation (6-30) is the final form for the general expression 
determining the storage changes to be made in order to lower 
costs, following the path of steepest descent. Again, a 
necessary condition for best economy is that (6-30) be equal 
to zero throughout the time interval considered. 

6.4 Transmission Loss in the Long Range Solution 

The determining equations thus developed do not take into 
account the effect of transmission loss on long range water 
allocation. This effect is present because some freedom exists 
in selecting the time when the water resources of a given plant 
should be used. It follows that this timing should be adjusted 
to reduce transmission loss if this does not materially increase 
the hydro deficiency at any time. If, for example, a hydro 
plant delivers its energy via a heavily loaded grid, the 
"transportation loss" of its own energy would be less if it 
were delivered during a week having lower average system load. 

Transmission loss can be added to the formulation devel¬ 
oped above. The additional complexity is appreciable, however, 
and the necessity for this addition is still in doubt.^ Formulae 

^•A study of the relative changes in operating costs due to in¬ 
clusion of transmission loss factors and changes in predicted 
stream flow has been undertaken by P. Johannesson in a companion 
M.S. thesis. 



- 62 - 

for the transmission loss in terms of weekly average genera¬ 
tions have been developed (Weiner and Whitbeck, 67) but re¬ 
quire further evaluation. 1 While this phase of the develop¬ 
ment is very incomplete, a formulation including transmission 
loss will be presented in order to indicate the order of mag¬ 
nitude of additional complexity incurred. 

Again, it is necessary to attempt simplifications of the 
formulation consistent with system operation. In the follow¬ 
ing, attention is focused on the load allocation among hydro 
plants, and an assumption is made that all sources for replace¬ 
ment of hydro deficiency can be characterized by a single cost 
. .2 

characteristic and loss coefficients for a single equivalent 
source. (The latter assumption will introduce errors depending 
on the distribution and size of non-hydro sources, but should 
be allowable for systems with relatively little steam). Adding 
transmission loss, equation (6-24) becomes: 

“ Ltt> + - £ *‘ (t) ( 6 - 31 ) 

and 

9&d _ 9?l _ 

?X ~ 9* L 9K ( 6 - 32 ) 

K*/ 

The loss formula developed by Weiner and Whitbeck (67) from 
eq. (2-1) is: 

■^A study of the accuracy of loss formulas has been undertaken 
by Mr. E. Greve in a companion M.S. thesis. 

2 

section 6.1 
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w»i n+t 




(6-33) 




where K(t) is a nearly constant factor taking into account 
the deviation of hourly loads from the weekly average. The 
A and B’s are constants depending on the transmission system, 
Separating out the terms pertaining to non-hydro sources: 


Pl = £ /L Sjf 2" + y 2 ®JW £<t+B<J4& + 

Jsl nt: I jfTJ 


L(i) K(t) 


(6-34) 


Using (6-34) in (6-31), one gets: 

ft[ ft] = *■«[ i + -#1 

L *•* J L M J (6-35) 

+ rr ft ft. 

n*»i 

This equation must be solved to determine gd(t). It is evi¬ 
dently more involved than (6-24) where transmission losses 
are omitted. It can be handled, however, and simplifies con¬ 


siderably when 


ft « i 2 ftd g, 


( 6 - 36 ) 


The expanded form of (6-32) is nsxt obtained. From (6-33): 


n *i 

Ml. =y 

'ft L 


Z B HJg 


Hence: 


-t-[ 

K=l /=! J j 


(6-37) 


( 6 - 36 ) 
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Using this in (6-32) one gets: 


iff 1- *!--£[ I 


3*. 


JX 


(6-39) 


In general, then: 

n ntt 

3g<(t*v t ) _ _ Y" r I pg«»Q 

^ 4 -^' ^ 

From the physics of the problem, it is knovm that: 

Ct+Vj) __ /-) K « 

P CO 6 *‘ 

^ & (t + T;k,Y __ Q K ^ b 
<? £ CO 

Therefore, using (6-41) and (6-40), and neglecting penalty 
functions, the two terms of the basic determining equation 
(6-3), are: 

PjSCt) 

as,ct) (6 - 42) 


jjjj4r (fc 0!) _ _ yfeffl) 

s 5. co 


i*< _ 


-f 2B^gyft) 




(6-40) 


(6-41) 


and: 


<j_r“ a$/hr(sAt+T4 

<** 1— 3 SM) 

0 r * 



(6-43) 

(t+Tb) 

a s,- ct) 


These are the same as those appearing in (6-25) except that 
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each % is multiplied by a correction factor. It follows, 
therefore, that equations (6-30) are also the determining 
equations for the case involving transmission loss, if the 
incremental cost, ^ , is replaced by: 

X‘(2Jnr,$ 5 ir(g,,&0 




n+i 


( 6 - 44 ) 


This correction plus the fact that gd must be evaluated using 
the more complex form (6-35), makes (6-30) applicable to a 
more general class of problems including transmission loss, 
but under the assumption that a single effective source of 
hydro-deficiency replacement can be formulated. 



Chapter VII 

SOLVING THE DETERMINING EQUATIONS 

7.1 Introduction 

The problem remaining is to solve a set of ordinary 
differential equations, the determining equations in ( 6 - 30 ), 
subject to boundary conditions at the beginning and end of 
the time interval to be optimized. The n um ber of variables 
and the complexity of the functions involved make the mecha¬ 
nization of the solution imperative. 

7.2 Digital vs. Analog Computers 

Both digital and analog computers are capable of solving 
this type of problem. 1 Digital computers have the advantages 
of a) higher accuracy obtainable by increasing the number of 
digits carried in the computations, b) versatility, in that a 
general purpose machine capable of performing other management 
and operational functions can be used, and c) size stability, 
for a single arithmetic element is used for all operations, 
and does not increase in size with growth of the system. As the 
system grows, some increase in digital storage capacity would 
be desirable in order to speed up the solution. However, the 
analog computer must add components corresponding to each 
system addition, thus increasing in size and multiplying cumu¬ 
lative errors. A further disadvantage of the analog type is 

Information on the use of the differential analyzer for solving 
the three plant problem will be included in the Master's Thesis 
by P. Johanneson. 


- 66 - 



- 67 - 

the inability to satisfy future boundery conditions for each 
plant except by successive trials. This is countered, of course, 
by the argument that the digital solution is inherently itera¬ 
tive. However, in the digital solution proposed, each iteration 
gives a useful solution that satisfies boundary conditions, and 
is an improvement, economically, over former steps. This type 
of solution lends itself particularly to cases where modified 
solutions are periodically called for corresponding to changes 
in system resources or loads. For these reasons, emphasis has 
been placed on the use of digital procedures for solving the 
determining equations. 

7.3 Philosophy of Approach 

Equation 6-30 in effect is an evaluation of the gradient 
of the performance surface^, telling in what direction the 
operating point should be moved or what changes should be 
made in storage elevations throughout the interval being opti¬ 
mized. The computer must evaluate this expression to modify 
operating curves and to improve economy by the path of steepest 
descent. 

The operator is to submit to the computer a set of drawdown 
curves representing a first approximation to good operating 
practice, together with stream flow and load predictions. The 
computer is then to automatically successively modify these 
proposed curves, giving with each modification an improved 
operation with reduced effective cost, while maintaining the 
initial boundary conditions. The movement of the operating 

^"Section 4.7 




point by the path of steepest descent involves the coordination 
of all storage changes at all plants throughout the time 
interval being optimized, so as to obtain the maximum possible 
reduction of the cost index for each increment of "integrated- 
storage-change-squared" given by: 

tl f ss- (if dt 

(-1 " c > 

However, since the curves operated on are relatively smooth, 
changes are made at a finite number of samples, say weekly 
samples, of the storage curves of each plant. 

Also certain approximations are used to evaluate (6-30) 
and finite rather than infinitesimal steps will be taken on 
the performance surface. Hence the path of steepest descent 
will not be followed perfectly. The essential feature is that 
each modification to system operation gives a more economic 
operation, while satisfying boundary conditions and avoiding 
or reducing the violation of operating limitations. 

7.4 Whirlwind I Computer 

The gradient method of solution involves repeating a. 
fixed sequence of calculations on successive sets of variables, 
namely stream flows, system load, and storage elevations at 
selected time intervals. This requires 1) a high speed computer 
which can execute each sequence of calculations in a relatively 
short time and 2) efficient data handling equipment which will 
rapidly and conveniently supply the computer with the successive 
sets of variables to be operated on. 
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The Whirlwind I digital computer has been made available 

for the economy loading study, and hence the demonstration of 
the gradient method is described in WWI terminology. However, 

any high speed digital computer is capable of performing this 
optimization with more or less speed and convenience, depending 
on its particular features. WWI computes at unusually high 
speeds, performing an average of ten thousand operations such 
as addition, subtraction, multiplication and division in one 
second. However, it is at present somewhat deficient in data 
handling or input-output devices. Hence the time spent by WWI 
in carrying out the sequences of calculations is negligible 
in comparison with the time spent in feeding data into the 
computer and obtaining printed results. Equally fast results 
could be obtained on a slower computer with faster in-out 
equipment. Before the delayed print out by magnetic tape was 
made available, each step of the iteration required about 15 
minutes of WWI time to execute. Additional time was required 
to transfer the results from punched paper tape to typewritten 
sheets. With the magnetic tape output, each iteration will 
consume less than five minutes of computer time. Again, most 
of this time is needed for data handling. 

A photograph of the control panel of WWI is shown in 
figure 7.9. Here the progress of the program is indicated and 
the location of component or program failures are presented. 

A photograph of the magnetic tape data handling equipment is 
shown in figure 7.14. 
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7.5 Computer Control Cycles 

Obtaining a solution on the WWI digital computer is 
divided in three parts: 

1) All the data that is to be used in the solution must 
be put in a form that is readily accessible to the computer. 

In this case, one must record on magnetic tape the values of 

the stream flows, system load, and original storage elevation 
curves, for each of fifty four weeks. 

2) Every step of the data handling and computational 
procedures must be coded on paper tape and these instructions 
must be delivered to the electrostatic storage of the computer. 
The computer then carries through the entire solution auto¬ 
matically and prints the desired results in final or inter¬ 
mediate form. Final form in the hydro problem is now typewritten 
results. However, present internal storage limitations reo^ire 
that an intermediate form, either punched paper tape or magnetic 
tape, be used. 

3) Results on an intermediate medium must be converted 
and transferred to final form. 

a) Magnetic Tape Layout 

At the top of figure 7.2 is shown a sketch of the magnetic 
tape with fifty-four blocks of data laid out, one for each 
week of the year being optimized, plus end weeks. Each such 
block contains twenty numbers, of which some are system re¬ 
sources and load for that week and some are the results of 
intermediate computations for that week. Initially only the 
former are recorded, with space left for the latter. The 
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former include: F^(t), F-^t + iTi ), F 2 (t), F 2 (t 4TT* ), L(t), 

L(t + *i»), y ^(t) andy 2 (t). When read by the computer, these 
are first transferred to electrostatic storage registers 
designated Sblthrough 22bl in Appendix F.lS.^" When a block is 
re-recorded by the computer during the optimization, the 
contents of registers bl through 3&bl are transferred to tape, 
thus laying out the results of intermediate computations as 
well as the original data. Each block is preceded by a block 
marker which enables the computer to locate each new block. 

Searching or skipping is achieved by counting block markers. 
Each block occupies approximately five inches of magnetic tape. 

b) Main Flow Chart 

A flow chart of the routine followed by the computer in 

carrying out the solution shown in figure 7.1. The order-by- 

order cycle-control program which executes this is given in 

Appendix F.2. As shown, two passes of the magnetic tape, over 

the fifty-four blocks of information, are needed to obtain 

2 

one step of the iteration process, or one ’’run”. The number 

of such runs to be performed by the computer is preselected, 
and the "run counter” in figure 7.1 keeps the computer on the 
job until the called for number of runs is completed. The "pass 
counter" keeps track of whether the magnetic tape is being 
passed for the first or second time in each run. The ordinate 
counter counts the weeks that are processed during each pass. 

^Each number on the tape occupies two 16 binary digit registers 
in electrostatic storage. 

2 

In the following the terms iteration, step, and run are 
used interchangeably. 





RESET PASS COUNTER 
AND ORDINATE COUNTER 


READ 1 GROUP 

SAVE PARTS IN "b" REGISTERS 
INSERT 2At IN 60bl 


SKIP 1 GROUP FORWARD 
READ 1 GROUP 

SAVE PARTS IN "f" REGISTERS 
SKIP BACK 2 GROUPS 
READ 1 GROUP 

IS PASS COUNTER POSITIVE? 


ES 


PRELIM. FINAL 
COM PUT. COM PUT. 




ES 


IS RUN COUNTER POS.? 


SKIP 1 GROUP FORWARD 
READ IN REPRINT 
INSTRUCTIONS 
REPRINT - END 


SKIP BACK 53 BLOCKS 



IS PASS COUNTER POS.? 


YES 


(D 


IS ORDINATE COUNTER > 1? 


YES 



EVALUATE EFFECTIVE COST 
ADJUST a 
INSERT At IN 60bl 
READ 1 GROUP 

SAVE PARTS IN "f" REGISTERS 


STEP PASS COUNTER 
SKIP BACK 53 BLOCKS 
RESET ORDINATE COUNTER 


FIGURE 7.1 FLOW DIAGRAM OF COMPUTING PROCEDURE 




















-BLOCK MARKER AND A BLOCK OF TWENTY 

- NUMBERS RECORDED ON MAGNETIC TAPE 

RD READ ONE BLOCK OF DATA FROM MAGNETIC TAPE AND 
TRANSFER TO ELECTROSTATIC 3TORAGE 

S SKIP OVER RECORDED DATA 

RC RECORD ONE BLOCK OF DATA ON MAGNETIC TAPE 
FIGURE 7.2 DATA HANDLING ROUTINE USING MAGNETIC TAPE 








- 74 - 


c) Data Handling Routine 

Although different computations are made and different 
quantities are re-recorded on the magnetic tape during the 
two passes, the same cycle of tape operations occurs in both 
passes, as indicated in figure 7.2. The starting and end 
points of the two passes are also different, the second pass 
omitting the end ordinates. The tape cycle of figure 7.2 
serves to a) bring to the computer data needed for the com¬ 
putation of corrections at one week, say week n, and b) re¬ 
record the data of block n including either intermediate 
results or modified storage elevations, depending on whether 
a first or second pass is being executed. 

Data needed for computations at block n include some 
stored at adjacent blocks, n+1 and n-1, which are used to 
obtain derivatives with respect to time. Data so used on the 
first pass are the storage elevations; and data so used on 
the second pass are computed values of the partials: 

jgkrCUTh) 

3 S; (t) 

The routine followed by the computer in calling for and 
temporarily storing magnetic tape data, as outlined in figures 
7.1 and 7.2, proceeds as follows: 

1) If at the beginning of a pass, read in the first block 
of that pass (block 0 for first pass; block 1 for second pass) 
and save portions of this block (needed for time derivatives) 
in electrostatic storage designated as ,? b" registers. If not 
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at the beginning of a pass, after computing for block n-1, 
save the corresponding portions of block n-1 in b registers. 
This stores certain data concerning the week prior to the one 
where corrections are next to be made. 

2) If not at the 53rd block: a)skip forward over the 
next block, n, and read forward the one following, n+1. Save 
the corresponding time-derivative portions of this block in 
the ”f” registers, b) skip backwards over two blocks, n+1 
and n; and read forward block n. Thus complete data for 
block n and time-derivative data from blocks n+1 and n-1 are 
then available in electrostatic storage. 

If at the 53rd block: read forward block n and trans¬ 
fer time-derivative portions to "f” register. Thus complete 
data for block n and certain data from blocks n and n-1 for 
time derivatives are available in electrostatic storage. 

3) After computations are made, skip back over block n 
and re-record forward block n, containing both the original 
stream flow and system load data for week n and also either 
the intermediate results of pass one or the storage elevation 
corrections of pass two. 

The above-described data handling cycle is repeated for 
each week in the interval being optimized. Operations progress 
from week to week as indicated in figure 7.2, in effect ’’level¬ 
ing-off the economic peaks” at each week in proportion to its 

magnitude. This recycling, indicated by the inner loop in 
figure 7.1, is achieved by a) adding one, after each week has 



- 76 - 

been processed, to a counter which initially has been set to 
minus fifty-one, b) recycling so long as this counter remains 
negative or zero. 

d) Examination of Cost Index 

After weeks one through fifty-two have been processed on 
the first pass, making the ordinate counter ♦!, the cost index 
is examined approximately by summing up the dollars per hour 
average costs for weeks one through fifty-two. (This cor¬ 
responds to a rectangle approximation for each weeks' area 
under the dollars per hour curve, so that the summation of 
weekly dollars per hour times 163 hours per week gives an 
approximation to the annual cost index. 

e) Adjustment of Size of Step 

At this time, the change in cost index achieved by the 
last iteration is also examined. If the change is negative, 
and is greater in magnitude than a preset amount, say $5000, 
the factor is left unchanged. If the change in cost index 
is negative, but is less in magnitude than, say $5000, the 
factor o< and hence the "size of step” is doubled. If the 
change in cost index is positive, indicating an overshoot of 
the performance surface minimum, the factor oC is halved. 

f) Recycling 

Following this procedure, the preliminary computations 
of the first pass are then completed for the 53rd block, as 
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indicated by the central return loop in figure 7.1. After 
this is done, the ordinate counter is again stepped (add one) 
making it *2. As shown in figure 7.1, the path is then clear 
to step the pass counter, making it *1, reset the ordinate 
counter to -50, and recycle via the outer loop to begin the 
second pass. 

The second pass inner-cycle in figure 7.1 is like the 
first, except that the positive pass counter ’’opens the 
channel” for the final computations instead of the preliminary 
computations and, if the run counter is still negative, re¬ 
cycles to begin another first pass after fifty-one ordinates 
have been processed. 

7.6 Computer Programs Used in the Demonstration of the 
Gradient Method 

The determining equations used to demonstrate the gradient 
method were the forerunners of equations (6-3), expressed in 
terms of elevations rather than storages.^ These are: 

-A- 

At 

where the summation must include the ith plant and all its 
downstream plants. With two hydro plants, dam 2 downstream 
from dam 1, one has therefore: 

■^Cypser, R. J. Engineering Memorandum No. 1, D.I.C. 7042, MIT, 
October, 1952. 
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a) Procedure 

The procedure used was to compute the dollars per hour 
cost for the given operation and also with small changes in 
each elevation and rate of change of elevation. That is, 
operating costs were computed at each weekly ordinate for the 
following conditions: 

TABLE 7-1 


1. 


i.(t) 

Hitt) 

Hitt) 

2. 

y,(t)+Au, 

if.Ot) 

i htt) 

Hitt) 

3. 

m 


iW) 

Hitt) 

4 . 

act) 

S.ct) 

%(t) + &Hi 

Hitt) 

5 . 

y,ct) 

S.(t) 

Hitt) 

Hitt)+ aHi 


The differences in operating costs thus obtained were used 
to approximate the derivatives in (7-2) and (7—3) * 

^Because of truncation and round-off errors, this method of 
evaluating the partial derivatives is inferior to that described 
in 6.2 and 6.3 leading to equation 6-30. 
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A list of the computations used to evaluate (7-2) and 
(7-3) at each ordinate is given in Appendix E in the following 
groups: 

a) Slopes of the drawdown curves 

b) Total plant discharge 

c) Net head 

d) Turbine discharge 

e) Hydro generation 

f) Hydro deficiency 

g) Effective hourly cost 

h) Components of (7-2) and (7-3) 

i) Storage elevation corrections 


At each ordinate, each of the quantities (a) through (g) is 
computed for each of the conditions in table 7-1. Used in 
these computations are expressions for Hir(Qi) , ii. ? ])<•) , 

( ^»i) I 8t'max: * and $/ hr ^d^» shown in figures 5.2, 
5.3, 5.4, 5.5, 5.6, 5.7, 5.10, 5.11, and 5.12. Also used were 
incremental storage curves given by the slopes of the storage 
curves shown in figures 5.S and 5.9. The implementation of 
the computations listed in Appendix E was carried out using 
the computer programs listed order by order in Appendix F. 


b) Subroutines 

Groups of orders that are used repeatedly are catalogued 
in Appendices F.3 through F.16. The principal groups are sub¬ 
routines p and k, the preliminary computations and final com¬ 
putations, used in the first and second passes, respectively, 
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as shown in figure 7.1. These in turn call upon the other 
subroutines listed, some being called upon repeatedly for 
different arguments, as described above. In the following, 
only a brief description of each subroutine is presented, 
with complete information relegated to the appendices. 

(F.3) Subroutine f compiles effective hourly cost from the 
relation: 

■>/j|r ~t/hr(&o)+m % g M +n s g/ -f 51 +rp ih 

H Hi 

(F.4) Subroutine g is used to limit the size of elevation 

corrections when penalty functions or spillage occur. 
The necessity for this limiter has not been definitely 
established, but it was included as a safety device to 
damp out possible overshoots and oscillations. These 
may be caused by larger corrections at certain weeks 
where operating-limitations are violated or spillage 
occurs, coinciding with an abrupt steepness in the 
performance surface in certain directions. The limiter 
characteristic used is sketched in figure 7.3. 


used 



FIGURE 7.3 CHARACTERISTIC OF ELEVATION CHANGE LIMITER 




The limit |£!/ m { is of course adjustable, and such 
limiting was not used unless penalty costs or spillage 
existed, at the week being corrected, after one of the 
last two iterations. This was kept track of in coun¬ 
ters occupying registers bl through 6bl as shown in 
F. Id and subroutines m and n in F.7 and F.d. 


(F.5) Subroutine, h computes incremental storages by the 
relations: 


(F.6) 


$, = + m, s i/ ls + n , s cf'* 

= s l0 +* mu f/ 2S 

Subroutine k contains final computations for second 
pass, involving evaluation of: 


A 3t/hd$ . si 3t/krLt*T, z ) j jj/MQ 

At j £ (t) ’ At At 


Su(i) 


(imf 


E.O); 
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(F.7, F.d) Subroutines m and n evaluate 

2; = t: [ 1- ~ a.T -D.- 

and sum penalty costs. 

(F.9) Subroutine p contains preliminary computations for 
first pass, involving evaluation of: 

J$/hr(i) J$/hr(i+7ii) 2$/hr W 

MU) net) ; i CO 



(F.10) 


}$jkr (t+ ^7z) PJ|4r (■£) 2 $/hr (i) 

ai,(t) z y* a) 

Subroutine d directs the computer to read the next 
block of magnetic tape data into electrostatic 
storage registers bl through 3#bl. 

(F.ll) Subroutine s directs the computer to move the mag¬ 
netic tape backwards over N block markers, when the 
number N has been previously placed in the accumula¬ 
tor (computer arithmetic element). 

(F.12) Subroutine w directs the computer to delay further 

operations for a preselected time (3-20 milliseconds) 
to allow the magnetic tape to come to a complete stop. 

(F.13) Subroutine t transfers data needed for time deriva¬ 
tives from registers 20bl through 29 bl to either the 
"f" registers, 40bl through 49bl, or to the "b" 
registers, 50 bl through 59bl, depending on whether 
40bl or 50bl has been previously inserted in t3. 

(F.14) Subroutine r serves to print on magnetic tape inter¬ 
mediate and final results which later can be trans¬ 
ferred to final form as typewritten results. The 
routine records one nine decimal digit number, with 
sign and exponent, at a time, taking this information 
from the "multipl® register accumulator" registers 
1050, 1051 and 1052. Hence each result appears on 
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magnetic tape as three ’’numbers" preceded by two 
block markers to facilitate later transfer. 


(F.15) Subroutine z evaluates the storage elevation penalty 
cost by the expression 

P* = w » % 


where t/ x is zero unless a violation occurs and is 
equal to the amount of violation otherwise, i.e., 


Hx 


y - at 



The slope of the penalty cost, fixed by the coefficient 
Wy , is made different for costs applied to the two 
dams by previously inserting either or 

in 50bl. 


(F.16) Subroutine u evaluates the penalty costs for exceed¬ 
ing either a computed maximum allowable generation or 
a minimum generation taken as zero. The maximum 
allowable generation, as a function of net head or 
machine ratings, whichever is lower, is precomputed 
in subroutine m or n. The penalty cost is then given 

by . 

p g = m t e. 


where g x is zero unless a violation occurs, i.e., 



f g-e, 




0-6 


\ >0 
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It is noted that an equivalent penalty cost as a 
function of turbine discharge would be preferable 

in this case, to facilitate use of lower limits on 
discharge for navigation and fish life. 

7.7 Use of the Gradient Method on the Three Plant System 

The procedures described in 7.5 and 7.6 were applied to 
the three plant system shown schematically in figure 4.5, 
using the computations listed in Appendix E, the character¬ 
istics shown in figures 5.2 through 5.12, and the system 
parameters in Appendix F.l. Operation during one complete 
year, starting with the last week in July, was considered, 
using the stream flows and system load shown in figure 7.12. 

The computer was provided with a proposed set of draw¬ 
down curves, sampled each week, shown in figure 7.13. By the 
programs discussed above, the computer was directed to eval¬ 
uate the effective cost of the proposed operation, and to make 
small modifications in both storage elevation curves throughout 
the year in order to reduce the effective operating cost, 
coordinating all such changes so as to move the "operating 
point" to one of lower cost by the path of steepest descent. 
Because the proposed elevation curves would, if used, involve 
the violation of project limitations during the spillage sea¬ 
son, the reduction of effective cost had to include the reduc¬ 
tion of the cost of replacement energy and the reduction of 

limitation-violation penalty costs. 

The successive reductions of effective operating cost, in 
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each of five such small steps of the gradient method, are 
tabulated in table 7-1. Operating costs for both the first 
forty weeks, during which no penalty-costs occurred, and the 
full fifty-two weeks are tabulated. The reductions of re¬ 
placement energy costs in the first forty weeks, and the 
reduction of limitation-violations penalty costs, which 
predominate in the latter twelve weeks, are substantial at 
each step. Due to these five modifications, the reduction 
in fuel costs for the first forty weeks is $167,9&£. or 
1.45$ of the original fuel cost for this period. The yearly 
reduction in effective cost, consisting mostly of reduction 
of limitation violation penalty costs in the latter period, 
is $1,0&5,397 or £.36$ of the initial effective cost. This 

higher percentage is indicative of the fact that modifica¬ 
tions are made throughout the year in proportion to the 
savings that a given change will net, so that limitation 
violations and peaks in the fuel costs will be more quickly 
reduced. 

Each successive step tabulated in table 7-1 gives a 
substantial saving without indicating the approach of a 

point of diminishing returns or a minimum in the performance 
surface. Further information on the nature of the surface 
being traversed can be obtained by examining the surface 
gradients tabulated in table 7-II. The gradient method has 

been set up to make storage changes at each ordinate in pro¬ 
portion to the value of the component of the surface gradient 



TABLE 7-1 



SUCCESSIVE REDUCTIONS 

OF EFFECTIVE 

COST 



IN 

EACH OF FIVE STEPS OF 

THE GRADIENT 

METHOD 



Original Cost 

First 40 Weeks 
$11,522,512. 

Saving 

52 Weeks 
$12,971,719. 

Saving 




$33,649. 


$172,981. 


First Modified Cost 

$11,433,363. 

$35,901. 

$12,798,733. 

$396,556. 


Second Modified Cost 

#11,452,962. 

$34,354. 

$12,402,182. 

$ 88,280. 

1 

Third Modified Cost 

#11,413,60S. 

$32,376. 

#12,313,902. 

$202,761. 

0^ 

ON 

1 

Fourth Modified Cost 

#11,336,232. 

$31,708. 

$12,111,141. 

#224,819. 


Fifth Modified Cost 

#11,354,524. 


$ 11 , 886 , 322 . 



Total Cost Reduction 


$167,938 
or 1.45% 


#1,085,397. 
or 8 . 36 % 
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TABLE 7-II 





GIVING RELATIVE 

MAGNITUDES 

OF THE SURFACE 

GRADIENTS 



IN VARIOUS DIRECTIONS 


Step 

£ss: 

fjs. 

52 

z*s: 

*»38 

52 

zts: 

MS* 



Dam 

1 


1 

1.5766 

0.6331 

394.56 

396.73 

2 

1.0439 

0.3430 

275.60 

277.69 

3 

0.S494 

0.9179 

133.06 

169.63 

4 

0.5969 

1.0165 

154.26 

155.67 

5 

0.5396 

1.1116 

336.02 

342.09 

6 

0.5020 

1.2319 

393.25 

394.96 



Dam 

2 


1 

5.791 

1.403 

1572.5 

1579.7 

2 

5.361 

1.623 

1123.6 

1131.0 

3 

4.636 

1.652 

1064.5 

1070.1 

4 

4.692 

2.262 

1266.3 

1273.2 

5 

4.655 

2.551 

1396.1 

1405.5 

6 

4.740 

2.776 

1672.3 

1630.3 



Both Dams 


1 

7.366 

2.096 

1967.1 

1976.4 

2 

6.424 

2.671 

1399.6 

1403.7 

3 

5.637 

2.570 

1252.6 

1259.9 

4 

5.239 

3.279 

1420.6 

1429.1 

5 

5.395 

3.663 

1766.1 

1747.6 

6 

5.242 

4.003 

2066.1 

2075.3 



- S3 - 


along the corresponding axis of the space having weekly storage 
as coordinates. Therefore the computed values of storage change 
will be used as a measure of the components of the surface 
gradient. Since the space is orthogonal, when motion is con¬ 
sidered having components only along certain axis, i = 1 , 2, 

...n, the surface gradient in this restricted direction is 
given by 

I V$\ = K 51 S 

K - ( 

Accordingly, at each step described above, a measure is given 
in table 7-H of the surface gradients in several directions, 
corresponding to motion along the axis of the first twenty-nine 
weeks, the next seven weeks, the following fifteen weeks, and 
the entire year, for the ordinates of dam 1 alone, dam 2 alone, 
and both dams. It is seen from table 7-H that with each step 
the surface gets flatter in the directions of the first twenty- 
nine ordinates, indicating diminishing returns from further 
such changes. However, while the surface levels in the first 
four steps in the direction of the last fifteen ordinates, it 
becomes steeper again in the following two steps, indicating 
increased savings to be made by further corrections at these 
weeks. Also, in the intermediate seven week period toward the 
end of the drawdown season, the slope quite steadily increases, 
promising further savings. 

7.3 Nature of the Modifications 

In the following, data is presented which indicates how 
the optimization works, in affecting weekly costs and storage 



elevations. The method involves a gradual reshaping of the 
cost curve for the entire year, making the larger modifications 
in storage where changes are most beneficial to cost reduction. 
The result is that the method inherently attacks peaks in the 
cost curve, cutting down the tops and eating into the sides. 

The numerical methods used to evaluate derivatives are imper¬ 
fect, particularly when rapid fluctuations are present, so 
that with rapid fluctuations in the cost curve one observes 
a levelling of the peaks and an inadvertent secondary filling 
of the valleys. 

The hourly effective cost of operation with the original 
storage elevation curves is shown in figure 7.4. The peaks 
during the last twelve weeks are attributable to limitation 
violations. The time scale corresponds to that given in 
figure 7.12. The smoothing action of the gradient method is 
clearly demonstrated in figure 7.5b, showing the modifications 
made by five small steps of the gradient method during weeks 
eleven through thirty-one. Here the relatively large high 
frequency components are cut down and smoothed out. The cor¬ 
responding modifications to the cost curve in the first nine 
weeks and weeks thirty-four through thirty-nine are shown in 
figures 7.5a and 7.5c. These sections are relatively smooth 
and the modifications demonstrate an ’’eating away” of the 
sides of very smooth cost curves. 

A closer look at weeks twenty-three through twenty-nine 
is given in figure 7.6, where the original cost curve and all 
five modifications are plotted. The consistent, stepwise 
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reduction of the two sharp peaks to a lower, broader peak is 
clear. The same procedure that cut down each peak would 
thereafter reduce the remaining single hump in like manner. 

The nature of the terms in the Euler expression for dam 1 
during this period is illustrated in figure 7.$. Evident 
are the relatively small, well behaved nature of 3 S/ht(i)/ 

9 Hi (*) and the larger, oscillatory nature of the partials 
with respect to slope. The latter largely determine the sign 
of the elevation corrections. The term 3$/kr(t+T„)/ 

s y. U) was found to be completely negligible and was not 
plotted. The consistent stepwise revisions of the storage 
elevation curve of dam 1 during this same period are shown on 
an expanded scale in figure 7.7. The modifications at dam 2 
were very similar. 

The reductions of effective cost during the last twelve 
weeks was due primarily to the reduction of limitation viola¬ 
tions in the form of overloaded machines at times of high flow. 
This condition was caused by a) inadequate drawdown with early 
refill, and b) inadequate height above crest to spill high 
flows after dam was refilled. The successive modifications 
served to simultaneously, stepwise, relieve both conditions, 
as shown in the storage elevation curve for this period in 
figure 7.11. The drawdown is increased, and refill delayed. 

So long as the dam is full at a time of high flow, increased 
spillage is needed to reduce exceeding the maximum turbine 
discharge. Therefore, elevation is increased temporarily 
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during that period. Continuation of the process of delaying 
refill would cut into this elevation peak and gradually 

eliminate spillage. This occurrance is indicative of the man¬ 
ner in which "local” conditions, in the vicinity of the week 
being optimized, determine the change to be made, and suc¬ 
cessive iterations stepwise "propagate" relatively remote 
effects into a local area. This is basic to the method. 

Changes made in the storage elevation of dam 1 are shown 
in figures 7.10a, 7.10b, and 7.10c. The elevation reductions 
in the first ten weeks shown in figure 7.10a are due largely 
to the fact that during this period the dam was nearly full, 
and the small explorations, Atf, , caused spillage. Hence 
elevation was reduced, providing a margin of safety of Ay, ft 

below crest. Successive corrections become smaller and smoother, 
indicating the approach of an equilibrium condition in this 
region. Elevation corrections during weeks twelve through 
twenty-five, shown in figure 7.10b, are characterized by the 
oscillations described above superimposed on a trend towards 
reducing the amount of drawdown in this region. The amplitudes 
of the oscillations become successively smaller with each itera¬ 
tion, corresponding to the leveling of the peaks in the cost 
curve, described above. In the later portion of the drawdown 
season, shown in figure 7.10c, the trend in elevation corrections 
is such as to increase drawdown by increasing amounts, indicating 
the propagation of effects into this region which make increased 
drawdown more advantageous. 
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Numerical tabulations complementing the curves presented 
above are given in tables 7-III, 7-IV, 7-V and 7-VI. 

A study of truncation and round-off errors arising in 
numerical evaluation of time derivatives in the Euler expres¬ 
sions is presented in Appendix G. The errors obtained using 
the simplest Lagrange formula are significant but are not a 
serious problem. The magnitude of errors can be limited, in 
a given application, by appropriate selection of sampling 
interval and approximation formula. The method tends to can¬ 
cel rather than accumulate such errors since each iteration 
evaluates anew the path of steepest descent. Furthermore, 
as the solution progresses, the smoothing action inherent in 
the method reduces truncation errors. 

7.9 Future Work 

It is hoped that the study presented in this treatise 
will be a stepping stone towards efficient system evaluation 
and automatized system operation. Much, however, remains to 
be done. In the following, some of the more important areas 
requiring further study are listed. 

a) The investigation of the feasibility of using exist¬ 
ing commercially available digital computers to implement the 
evaluation of large systems by the gradient method. 

b) The introduction of probability distributions of 
stream flows in place of sets of expected stream flows, to 

organize operational optimization on the basis of "expected 


costs" 
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c) The evaluation of the importance of transmission 
loss, and the importance and feasibility of obtaining accurate 
predictions of stream flows, in the long range optimization. 

d) The use of the derived determining equations to study 
the relative value of changing certain plant parameters or 
capacities. 

e) Further development of methods of handling the short 
range optimization, coordinating this with the long range 

problem to obtain integrated procedures embodying both. 

f) Further study of the accuracy and utility of existing 

transmission loss formulae when basic assumptions made in 
their derivation are violated. 

g) Development of procedures or mechanisms to semi- 
automatically carry out the load dispatching of a system, 
based on long and short range economy optimization. 
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TABLE 7-III 

CHANGES IN EFFECTIVE HOURLY 

COSTS IN 

FIVE 

SMALL STEPS OF THE 

GRADIENT 

METHOD 

Week 

$/hr 

Week 

$/hr 

1 

4.3147 

27 

5.0206 

2 

-.2327 

23 

-34.6326 

3 

-21.3353 

29 

-31.0994 

4 

-19.5570 

30 

20.3609 

5 

-9.6513 

31 

-14.3231 

6 

4.3525 

32 

-4.4932 

7 

6.0224 

33 

13.1635 

3 

-16.1041 

34 

-2.9234 

9 

-3.6550 

35 

-31.6775 

10 

10.1499 

36 

-111.0235 

11 

17.3525 

37 

-242.4290 

12 

26.1250 

33 

-140.9379 

13 

-.4127 

39 

-362.7776 

14 

-33.2300 

40 

-43.2316 

15 

5.2933 

41 

-5.3317 

16 

62.5395 

42 

-192.3019 

17 

23.5763 

43 

-304.2360 

18 

-9.4204 

44 

-20.9225 

19 

39.9395 

45 

-330.4501 

20 

31.7313 

46 

-1160.3759 

21 

-53.2634 

47 

-1533.2603 

22 

-9.6960 

43 

-1247.2256 

23 

35.5472 

49 

-974.7123 

24 

-64.6423 

50 

-3.3722 

25 

-67.0715 

51 

-709.4002 

26 

33.3439 

52 

-43.3397 
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TABLE 7-IV 

AVERAGE CHANGES IN EFFECTIVE COSTS 
AFTER FIVE SMALL STEPS OF THE GRADIENT METHOD 


Weeks 


1 - 

5 

-$ 7,322.432 | 



6 - 

10 

$ 212.637 ' 

' -# 

4,934.27 

11 - 

15 

$ 2,625.521 ^ 



16 - 

20 

# 25,732.406 " 



21 - 

25 

-# 26,733.16$ 

-# 

9,611.43 

26 - 

30 

-# 3,660.669 . 



31 - 

35 

-# 6,347.462 ' 



36 - 

40 

-#151,275.533 ' 

f -#311,732.43 

41 - 

45 

-#153,609.490 _ 

1 


46 - 

50 

-#334,367.062 " 

I 


51 - 

52 

-#126,460.303 

( -#961,327.36 
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TABLE 7-V 

STORAGE ELEVATION CHANGES AT DAM 1 
IN SIX SMALL STEPS OF THE GRADIENT METHOD 


Week 


2 

-.0460 

3 

-.0207 

4 

-.0783 

5 

-.0800 

6 

-.0835 

7 

-.0535 

a 

-.0724 

9 

-.0808 

10 

-.0774 

11 

-.0565 

12 

-.0175 

13 

.0148 

14 

-.0148 

15 

-.0776 

16 

-.0037 

17 

.0994 

18 

.0949 

19 

.0732 

20 

.1842 

21 

.1566 

22 

.0453 

23 

.1334 

24 

.1175 

25 

-.0341 

26 

-.0337 


Week 


27 

.0705 

28 

-.0210 

29 

-.1499 

30 

-.1031 

31 

-.0953 

32 

-.1400 

33 

-.1073 

34 

-.1069 

35 

-.1162 

36 

-.1973 

37 

-.4288 

33 

-.9503 

39 

-.8834 

40 

-2.2818 

41 

-1.1428 

42 

-1.9211 

43 

-1.9867 

44 

-1.7774 

45 

.2682 

46 

1.4381 

47 

1.4621 

48 

1.9440 

49 

3.3683 

50 

3.4447 

51 

.7286 

52 

-.2681 



TABLE 7-VI 

STORAGE ELEVATION CHANGES AT DAM 2 
IN SIX SMALL STEPS OF THE GRADIENT METHOD 


Week 


2 

.0447 

3 

.0360 

4 

.0246 

5 

.0169 

6 

.0049 

7 

.0092 

6 

.0060 

9 

.0012 

10 

.0016 

11 

.0002 

12 

.0002 

13 

-.0013 

14 

-.0030 

15 

-.0026 

16 

.0003 

17 

.0042 

id 

.0039 

19 

« 

o 

o 

\_n 

-0 

20 

.0110 

21 

.0113 

22 

.0053 

23 

.0064 

24 

.0061 

25 

.0016 

26 

.0001 


Week 


27 

.0015 

26 

-.0009 

29 

-.0062 

30 

-.0047 

31 

-.0057 

32 

-.0063 

33 

-.0067 

34 

-.0069 

35 

-.0106 

36 

-.0220 

37 

-.0352 

36 

-.0676 

39 

-.0601 

40 

-.0904 

41 

-.0424 

42 

.1536 

43 

.0050 

44 

-.1997 

45 

-.2475 

46 

-.1037 

47 

.3002 

46 

.5310 

49 

.5656 

50 

.3932 

51 

-.2675 

52 

-.2661 
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and received his M.S. in Electrical Engineering in June, 1947. 
The summer of 1947 was spent at the Research Laboratory of the 
General Electric Company, Schenectady, N.Y. 
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Appendix A 

Path of Steepest Descent With Certain Coordinates Fixed 


In the following it is demonstrated that even though the 
coordinates corresponding to pertinent characteristics of the 
actuating inputs are fixed, maximum reduction of cost per 
length of step taken is obtained by making each change in the 
adjustable parameters equal (or proportional) to the correspond¬ 
ing partial derivatives of cost with respect to adjustable 
parameter change. 

Consider the space having the relevant characteristics of 
the actuating inputs as well as the adjustable parameters as 
coordinates. As the actuating inputs change,the operating 
point moves over the surface which is stationary. At each 
operating point however, a unique gradient exists on this sur¬ 
face whose direction determines the desirable change in the 
input variables. Of the input variables, only the adjustable 
parameters are controllable, so only the components of the 
surface gradient along these axes are of interest in the con¬ 
trol problem. The gradient at any point on the surface may be 
expressed as: 



terms dependent on 


where are unit vectors along the adjustable parameter 
axes, x^, are the adjustable parameters, and y^ are the rele- 
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vant characteristics of the actuating inputs. The gradient $ 
is a vector which points along the direction of maximum slope 
of the surface at the point where it is evaluated. This 
equation is rewritten: 


grad $ = IK* + b 


The vector b includes all terms not in the above summation and 
has no components along the adjustable parameter coordinates. 

An incremental change in operating point is equal to the cor¬ 
responding incremental change in the radius vector from the 
origin to the surface. Hence, the performance change correspond¬ 
ing to an incremental change in operating point is given by 
d$ * grad $ . dR 

The slope of the surface in the direction of the incremental 
change in operating point is given by 
d$/ds s Grad $ . dR/ds 

where s is the scalar distance moved along the surface, Separating 

gradient $ as above: 

df - 9$ . dR/ds + ~b . dR/ds. 
ds J X 

If the change in operating point is obtained by changes only 
in the adjustable parameters, dR is constrained to have com¬ 
ponents only along these axes. Since, b does not have any 
components along these axes, the latter term b.dR - 0, and 
dl/ds = . dR/ds (with constrained dl?). Since dR/ds is a 

unit vector, the dot product is a maximum when the two vectors 
are parallel. It follows, therefore, that under this constraint 
the maximum possible slope will be followed if the direction 



- 120 - 


of motion is parallel to the vector, -j-f . The relative 
magnitude of the components in this vector give the relative 
improvements in performance to be gained by a change in each 
of the adjustable parameters. The algebraic sign of each 
component determines the direction in which it would be advan¬ 
tageous to adjust each parameter. Regardless of how the 
operating point moves, because of changes in the actuating 
inputs, and regardless of possible malfunctions in previous 
steps of the path-finding process, each evaluation of the new 
gradient gives a fresh start toward the minimum. 
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Appendix B 

First Proof of Determining Equations for Three Plant System 

This first proof is obtained by a geometrical interpreta¬ 
tion in a multidimensional, orthogonal, coordinate system. 

This approach enables one to visualize the optimizing procedure. 
Let there be: 

1 sampled ordinates of S-jJt) in the interval 
m sampled ordinates of S 2 (t) in the interval o N < t 4 T 
Equation (5-37) is then expressed as: 
j Jlim 

s$=y /s,(vE,(tut +y / 5, ct ; ) E,ctjAt + < (B _ 1 j 

i-.JU l 


where 


E,tt) 


A 96/hr( t) A 

»s,m -Jt is,(tj-Jt P 3 , (t;] < B - 2 > 


Ej (£•) e il lkM _ jL iMt (*1 
J Si (t[) dt 3 Si(tf) 


The number, €, can be made arbitrarily small by letting 1 and m 
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be finite but arbitrarily large. Therefore even neglecting 
6 in (B-l) t can still be found to any desired accuracy by 
increasing 1 and m. 

Although the cost index is also dependent on the physical 
and economic characteristics of the system and on the pre¬ 
dictions of stream flows and system load, only the storage 
curves S-^(t) and S 2 (t) are controllable and are here considered 
as variables in the optimization. For the case considered here 
where the cost index is made a function of a large but finite 
number (1 ♦ m) of storage ordinates, ordinary calculus gives us 


6 


A 



i = i 


n 

) $, ft) 


<r s, ft) 


}S t (to 


S s/t) 


(B-4) 


If 1 and m are large enough so that 6 may be neglected in (B-l) 
then the forms of (B-l) and(B-4) are identical. It follows, 
then, that: 


i$__ 

d %(t^) 


E*1 (t^) A t 


/ ^ i 4 4 (B-5) 


9 f 


E2( ti) A t 


4 i * J? i m (b-6) 


A vector space is now defined, containing: 

1 orthogonal axis S^('ti) i z 1, 2, ...% 

m orthogonal axis S 2 (ti) i r 1+1, 1+2, ...i*m 

Corresponding to each point in this space, there are, as its 
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projections, specific values of the S^(ti) and S 2 (ti), and 

hence there is a specific value of the system’s cost index. 

a 

Associated with each axis is a unit vector: V t '■ *»*»•• Jf+n* 

A small step taken in this space is expressed as the 
change in the radius vector: 


je+m 


<rR =^T S S, (t) $■ +Y_ sS,(i) 


(B-7) 




i=J+t 


Also, the gradient of the cost index is defined as: 

lltn 


Vft = 2 _ 


^ n 


777 


7 V/ 


4- \ 3 % A 


i 

i*J +« 


V- 


3S i Ct) ‘ 


(B-S) 


Therefore, since* the dot product of two vectors is the sum of 
the products of corresponding components, equation (B-l), 
modified by (B-5) and (B-6), may be rewritten as: 


S$ = SR ' V$ 


(B-9) 


SR is an incremental line segment in space, along which the 
"operating point" moves, and whose direction depends on the 
relative magnitudes and the signs of all the S and £’S 2 . 

The problem is the selection of the direction of SR , that is, 
the relative changes in all the ordinates, so as to obtain the 
greatest reduction of the cost index for a given size of step, 

I/ll. 
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Although the cost index is a single-valued function of 
the storage ordinates, there are many sets of ordinates having 
the same cost index. Hence, since the cost index is a con¬ 
tinuous function of the ordinates, through each point in 
space there passes a surface of constant cost index. 

With a given size of step, the magnitude of & $ will be 
a maximum when Tt is oriented perpendicular to a surface of 
constant cost index and parallel to the gradient of the cost 
index. This is seen from (B-9), for the magnitude of a dot- 
product of two vectors with fixed magnitudes is a maximum when 
they are parallel. To make S $ negative, & R should be in a 
direction opposite to ~n. This requires that the components 
of S R be equal to the corresponding components of 7 $ except 
for a common negative scale factor. That is, one needs: 



S S t ft). - - 


3 $ 

9 Sz(t;) 




(B-10) 




Where °< is a positive constant determining the size of step. 
Using (B-5) and (B—6) in (B-10), the desired changes are given 
by: 


SSA)*- 


MjrCtd . _ -V sMrCUti 

Jt ss,(t : ) clt 


& t 


(B-ll) 


i -J 
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ss t a)--* 


2t/hr(ti) 

3 S, Ct { ) 


a) 3 f/hr Ct) 
3 S 2 Ct) . 


At 


i = jH, ***•> 


(B-12} 


These equations are identical, to within a constant factor, 
to equations 5-40 and 5-41 and so complete the proof. 
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Appendix C 

Second Proof of the Determining Equations 
for the Three Plant System 


A second proof, which does not require shifting to 
finite sums, is given below. The objective, again, is to 
make $ $ of (5-37) as large, negative, as possible, with a 
limited magnitude of "integrated storage change squared". 

To proceed along the path of steepest descent, only 
changes in S-^(t) and S 2 (t) which definitely contribute to 
reducing cost must be made, and hence both integrals in 
(5-37) must be negative or zero. The problem is attacked 
by first making each integral of (5-37) as large, negative, 
as possible; then the relationship between the integrals 
is fixed. Equation (5-37) is expressed simply as: 


r 

sS,U)E,a)Jt 


+ 


T 

f s S, ft) El ft) 

a 


At 


(c-i) 


where E-j_ and E 2 are the Euler-type brackets in the first and 
second integrals, respectively, of (5-37). The measure of 
storage change at plant 1 is given by: 

T-Vz 

(S/t.f = [ sS.it) jt 


(C-2) 
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It is now assumed fixed, with a magnitude to be later specified. 
Schwartz’s inequality states that: 





(C-3) 


It follows that: 



Using (C-2) in (C-4), one gets: 


■ T-t. 


S S,(t)E,(t)<Jt 



(C-5) 


It follows, therefore, that the first integral of (C-l) must 
lie in the interval 




Hence, the most negative value possible for this integral is 
the lower bound of the integral (C-6), namely: 
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- U (c-7 

The storage changes needed to obtain this lower bound are 
given by: 


S 3 .«) = 


s/f, 


fTlmt 


£T, (t) 


(C-3) 


This is easily demonstrated. Using (C-3) in the first inte¬ 
gral of (C-l), one obtains (C-7) as follows: 


r-r. 


“<f $.tt) = -S4, 


T-1,% 2 

F Ctf dt 


•r-r,. 


E,(t) 3 dt 



(C-9) 


In like manner, the second integral of (C-l) can be bounded; 
and this lower bound can be obtained by making storage changes 


S S t «) = 


-6 ft. 


JJzMTt 


E,ft) 


(C-10) 


Where the measure of storage change at dam 2 is given by: 
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(s/?R =/ sstctfjt (c . n) 

Equations (C-&) and (C-10) are sufficient to make each of the 
integrals in (C—1) as large, negative, as possible, for speci¬ 
fied values of S and There remains to be demon¬ 

strated that the coefficients of E^t) and E 2 (t) in (C-S) and 
(C-10) should be equal if the sum of these integrals is tobe 
the largest, negative, possible under a constraint on the total 
storage change. The latter is defined as: 

sSzCtfdt 

(C-12) 

— <f 

Substituting (C-S) and (C-10) in (C—1) one gets: 

S$ = - S /T, J~j' E, df dt ~ Jf' (c—13) 

The objective is to determine the relationship between & 
and <f/f 2 which will make (C-13) stationary under the sub¬ 
sidiary condition that be a specified value, or that: 


&/? * = J fSXtfdt -t- J 
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= cf/f: 1 sfi* - sfi = o 

(C-14) 

_ l 

From the calculus , the necessary condition for 

this is that: 

— (£$L _j_ k ^ & -0 c-1** 

3 (3*;) a(sti) 

(C-15) 

where k is a constant, Lagrange’s multiplier. 

Using (C- 13 ) and (C-14) in (C-15), one gets: 


~ J f ~ E t (t)*dt -jr 2 K Sfi ~ 0 


f » 

(C-16) 

-jf E,(t?Jt + ZK f/f* = 0 


Solving these for and S ^ 2 : 


Srf, - i k f( E,(i) a dt 

(C-17) 

7k -/f *(*)'# 

(C-13) 

^Courant, 12 
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Substituting (C-17) in (C-S) and (C-lS) in (C-10) one finds 
that the coefficients of E-^(t) and E 2 (t) are equal to 



= — X 


(C-19) 


This establishes the relations for storage changes, so as to 
follow the path of steepest descent, as: 

s S,ct) = - <* E, (t) 

(C-20) 

s s g ct) - - * e 2 ct) 


thus completing the proof 
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Appendix D 

A Formulation of the Short Range Optimization 

Because cumulative variations of plant efficiencies 
during any one week can be neglected, the short range optimi¬ 
zation simplifies to an ordinary minimization in calculus* 

This is developed for the following case: 

1. Loadings of each hydro plant at any load below a 
predetermined capability limit are equally efficient. 

2. The cost of replacement of hydro deficiency is a 
known quadratic function for each source. 

3. Each hydro plant is constrained to deliver a speci¬ 
fied average weekly power. 

4. Transmission loss is appreciable. 

5. n hydro plants and m non-hydro sources comprise the 
system. 

The cost of supplying hydro deficiency from plant j is given 
by: 

$/kr (g d j) = f(o + 4 § Jt +-t fji V-m (d.!) 

The corresponding incremental cost is therefore: 

$/vw. hr (gjj) = = §1 + 


y- f-m 


(D-2) 
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The undesirability of exceeding a predetermined loading on 
each hydro plant is accounted for by using a penalty cost 
which increases abruptly when this limit is exceeded. This 
ficticious incremental cost is given by: 


y»«kr(g kl ) - ~ * + (i~j ** 

= + fu gfu 


(D-3) 


where a - o for g^and a » ( for g h . > g h , m 
The incremental costs for one source of deficiency power and 
one hydro plant are sketched in figure D.l. 



Figure D.l Sketch of incremental cost of hydro deficiency 
and incremental penalty cost for overloading hydro plant. 

If the penalty cost curve is made sufficiently steep (suf¬ 
ficiently large a) only small violations of the generation limit 
can occur. , 

Sj< A 

The specification of an average weekly generation' 1 at the 


hydro plants requires that: 
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/ (fifci - ^a) <Jt - ° ' * 


(D-4) 


The necessity of meeting the load requires that: 


L (V + p L (i)-t g,; ~ V g h , = O 


r> ’■>> 


(D-5) 


We seek to minimize the effective cost given by: 


O/ne*) +E fWsy 

"o 1 r' 




dt 


(D-6) 


subject to the constraints (D-4) and (D-5). These constraints 

are included in the minimization by adding terms to the inte - 
grand of (D-6).-*- We then seek to minimize 

-T 


I * / FCt) dt 


D-7) 


where 


FCt) = Z + Z 

i*< i=< 

~ f ft*" ft. A l _ ^ ft)[ l + p t "5 ft/ * 5! 


(D-d) 


^"Methods of Advanced Calculus, P. Franklin. McGraw Hill Book 
Co., 1944, p.439-442. 
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where the A; are constants and A is a function of time, 
all definite but unknown. 

The necessary conditions for a minimum of the integral of 
(D—7) are the Euler equations: 


3 F __ 9 F __ 


=r 0 


If. _ J£L _ 




o 


<■ * I,i •••" 




ID-9) 


Since F is not a function of or g^, this reduces to: 


n 

M,- 


= o 


; = i.2. 


Using (D-2) 


44- = o > = * 

(D-3) and (D-3), (D-10) becomes: 


(D-10) 


Pi, + fa 8k 



i = *,■■■* (D-ll) 

fji + -a($ [ 2 B/,hi ^ r ° 

/ - ! t 8, ■■■ ry\ 

These, together with (D-4) and (D-5) involve 2H+W+/ 
unknowns and an equal number of equations, which can be solved 
by either analog computer or numerical iteration. 
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APPENDIX E 
TABLE OF COMPUTATIONS 

USED IN THE DEMONSTRATION OF THE GRADIENT METHOD 
IN A THREE PLANT SYSTEM 

using the following approximations with small T^: 

h (t) = h (t + r i2> = ' V 

y 2 (t) = y 2 (t + t 12 ) =y 2 (t-r 12 ) 

^(^(t)) = S l^l^ +T 12^ = s l^ y l (t " r i2^ 
s 2 (y 2 (t)) = s 2 (y 2 (t + r i2 )) = s 2 (y 2 (t - r i2 )) 

a) Slopes: 

Cl y^Ct) = y x (t + At) - y x (t - At) / 2At 

c 4 y 2 (t) = y 2^ ~ y 2^ / 2^1 

b) Plant Discharge: 

C6 Q x (t) = F 1 (t) - ^(t) 

C7 Q ] _(t, y x + Ay x ) = F 1 (t) - S^t, y 1 + Zty^) 

C8 Q x (t, y ± + Ay x ) = F 1 (t) - S^t, f ± + Ay^ 

C9 Q,(t + r 12 ) = F x (t + r 12 ) - s x (t) 

CIO Q 2 (t) = F 2 (t) - S x (t) - S 2 (t) 

Cll Q 2 (t, y 2 + Ay 2 ) = F 2 (t) - S-^t) - S 2 (t, y 2 + Ay 2 ) 
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C12 

Q 2 (t, y 2 + Ay 2 ) = F 2 (t) - S 1 (t) - S 2 (t, y 2 + Ay 2 ) 


C13 

Q 2 (t + r i2> ? ? 2 (t + V " - S 2 (t) 


C14 

« 2 (t + r u , y x + asrj) = r 2 (t + r 12 ) - s 1 (t,y 1 + ^ 

- 

C15 

Q 2 (t + T 12’ h + 5 F 2 (t + T 'l2 ) " + 

- s 2 (t) 

NOTE: 

s ± (t) = s 1 (y i (t)) y^t) 



c) Net 

Head: 


C16 

^(t) = y 1 (t) - y 1T (C6) 


C17 

y 1 + /^ ± ) = y x (t) + ^ - y 1T (C7) 


CIS 

^(t, y 1 + Ay x ) = y x (t) - y 1T (C8) 


C19 

h l^ + r l2^ = y i^ + r i2^ “ y lT ^ 


C20 

b 2 (t) = y 2 (t) - y 2T (CIO) 


C21 

h 2 (t, y 2 + Ay 2 ) — y 2^^ ^2 ™* y 2T 


C22 

h 2 (t, y 2 + Ay 2 ) = ~ y 2T ^ C12 ^ 


C23 

h 2 (t + r 12 ) = y 2 (t + r 12 ) - y 2T (ci3) 


C24 

h 2 (t + r 12 , y x + ^ = y 2 (t + r 12 ) - y 2T 

(C14) 

C25 

h 2 (t + T 12’ y l + = y 2 (t + r i2 ) “ y 2T 

(C15) 


d) Turbine Discharge: 


026 D 1 (t) = C6 -i 1 (y 1 (t)) 

G27 D 1 (t, 7 1 + Ay^) = C7 - J 1 (y 1 (t) + Ay^) 
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C28 D 1 (t, f ± + Z^ 1 ) = C8 

C 29 D 1 (t + r 12 ) = eg -^ 1 (y 1 (t + r 12 ) 

C 30 D 2 (t) = cio -,/ 2 (y 2 (t)) 

C31 D 2 (t, y 2 + Ay 2 ) = Cll - + 

032 D 2 (t > ^2 + ^2 )= 012 * J 2 ^2 {t)) 

C33 R,( t + r 12 > =C13-7 2 (y 2 (t + r ls )) 

C34 D 2 (t + r 12 , y x + % x ) = CU -i 2 (y 2 (t + r 12 )) 

C35 D 2 (t + r 12 , y x + Ay x ) = C15 - J 2 (y 2 (t + r i2 )) 


e) Hydro Generation: 


C36 g-^t) = ~ (C16) • C26 

g-, 

C37 g x (t, y;L + ZV-l) = jj“ (C17) • C27 


038 


g i 

Sj_(t, y ± + = ~ (CIS) * C28 

g q 1 

C39 g x (t + r i2 ) = ^ (C19) * C29 


C40 g„(t) = ^ (C20) ♦ C30 


g 2 

Ca g 2 (t, y 2 + Z^ 2 ) = ^ (C21) • C31 
g 2 (t, y 2 + Z^r 2 ) * ^ (C22) • C32 


C42 

C43 


+ r !2> = T 2 < C23 > ' 033 

So 


C44 g 2 (t + r 12 , 7l + 4y x ) = ~ (C24) ' C34 
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C45 g 2 (t + ? 12 , f ± + /^ ± ) = (C25) ' C35 

f) Hydro Deficiency: 


C 46 

g d (t) = L(t) - C36 - C40 


C47 

g d (t, y x + ZJy x ) = L(t) - C 37 - C 40 


048 

g d (t, y^ + Ay x ) = L(t) - C38 _ C 40 


C49 

g d (t + r 12 ) = L(t + r 12 ) - C39 - C43 


C50 

g d (t, y 2 + Ay 2 ) = i(t) - 036 - C 41 


C51 

g d (t, y 2 + Ay 2 ) = L(t) - 036 - C 42 


052 

g d (t + r i 2 > y i + ^ 1 ) = L (t + * 12 ) - C39 

- C44 

C53 

g d (t + r i2' y l + ^ 1 ) = L ^ t + r i2^ ** 039 

- C45 


g) Effective Hourly Cost: 

C54 $/hr(t) = |/hr (C46) + p lg (036) + p^y-^t)) 

+ P 2g ( 040 ) 

C55 $/hr(t, y ± + Z^) = I/hr (C47) + p lg (C37) + P ly (y 1 (t) + 4y x ) 

+ P 2g (C40) + 

C56 l/hr(t, y x + /V x ) = $/hr (048) + p lg (038) + P^y^t)) 

+ p 2g (C40) + p 2y ( y2 (t)) 

$/hr(t + r 12 ) = I/hr (C49) + p lg (C39) + p^y^t +r 12 )) 

+ P 2g (0*3) + P2y(y 2 (t + T 12 )) 


C57 
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058 |/hr(t, y 2 + 4y 2 ) = I/hr (C50) + p lg (C 36 ) + P ly (y 1 (t)) 

+ P 2g < CW) + P2y<y 2 (t > + ^ 2 > 

C59 l/hr(t, y 2 + & 2 ) = I/hr (C51) + P lg (C 36 ) + P ly (y x (t)) 

+ p 2g (C42) + P 2 y (y 2 (t)) 

C 60 l/hr(t + r 12 , yi + /Ny x ) = |/hr (C52) + p lg (C39) + P 3 y (y 1 (t +r 12 )) 

+ p 2g (044) + P2y (y 2 (t +r 12 )) 
C61 |/hr(t +r 12 , y 2 + Ay x ) = |/hr (C53) + P lg (C39) + P^y^t + ^ 12 )) 

+ P 2 g (C45) + p 3jr ( y 2 (t +r i 2 » 

h) Components of the Euler Expressions: 


C 62 

1 ! 

£ 

.c$iL^-C$£ 


ay-j/t) 

^1 

C 63 

3 $/hr(t + r 

12 ) C 60 - C57 


3 y x (t) 

- AJ X 

C 64 

d |/hr(t) _ 

C $.6 


a j x (t) 


C 65 

3 |/hr(t +r 

12 ) 061 - C57 


a y x (t) 

Ay x 

C 66 

3 $/hr(t) _ 

3* 

0 

1 

CO 

in 

0 


3 y 2 (t) 

ay 2 

C 67 

3 &/hr(t) 

C59 - C54 


a y 2 (t) 

^2 
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dt 3 y , 1 (t) 2At 

d p $/hr(t +T 12 ) f 65 - b65 

<*t 3 y- L (t) ~ 2At 


C70 d 


d t 3 y 2 (t) 


f 67 - b6j 
2At 


i) Elevation Changes: 


C71 &y 1 (t) =-—-r [- C62 - C63 + C 68 + C 69 ] 

s 1 (y 1 (t)) 


C72 &y (t) =-=-- 

s 2 (7 2 (t » 


[- C66 + C 70 ] 
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Appendix F 

Programs Used to Implement the Gradient Method 

on Whirlwind I 
for the Three Planet System 
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F. l 

DOUBLE LENGTH 
PARAMETER STORAGE 


TAPE NO. 2260 P 1 
( 24 , 6 ) 

Storage 

Register 

32/bl, 


+.0 

DITTO THRU 


TEMPORARY STORAGE 


74 r/ 74 bi, 

+3336.30 

y 10 = y i mln 

76bi 

+2.S00 

8 i(yi 0 ) 

7 &bl 

+0.02599 

m ls 

gObl 

+O.OOOO5663 

n ls 

&2bl 

+3072.0 

y it(o 1 = 0) 

g 4 bl 

+1.5012 

m lt 

66bl 

+ 0.034039 

“ n lt 

2Sbl 

+20.0 

h 10 

9 Obi 

+256.00 

h ll 

92bl 

+ 0 . 064 S 6 

u 

94 b 1 

-3^.1075 

m lg 

96bl 

-2.8626 

~ n lg 

9 Sbl 

+ 128.7 

g lm(h = h lx ) 

lOObl 

+1.2299 

mig m 

102 bl 

~o.0005136 

~ n lg 

104 bl 

+3556.0 

y lc 

106bl 

+0.69265 

m lL 



-H4 


lOSbl 

+ 0.541125 

n lL 

110b! 

+300 

g lm(rated) 

112bl 

+0.300 

^1 

ll 4 bl 

+0.03 


ll6bl 

+ 2 SS 3.0 

y 20 = y 2 min 

USbl 

+60.0 

8 2 ( y 20 > 

120bl 

+0 

m 2s 

122bl 

+2697.OO 

y 2t(Q 2 =0) 

124 bl 

+ 0 . 56 S 75 

m 2t 

126bl 

+ 0.003594 

~ n 2t 

12Sbl 

+ 5.0 

0 

C\J 

A 

130bl 

+ 164.0 

h 22 

132bl 

+0.0760 

2 

134 bl 

-2.SS5 

m 2g 

136bl 

*• 0 . ISO 

"* n 2g 

13Sbl 

+120 

g 2m(h g = h 22 

l 40 bl 

+ 1.2299 

m 2gm 

l 42 bl 

-O.OOO5136 

~ n 2gm 

l 44 bl 

+ 2 S 93.0 

y 2c 

l 46 bl 

+0.27706 

m 2 L 



l 4 gbl 

+0.21645 

n 2L 

150bl 

+150 

^2m(rated) 

152 bl 

+0.1 

Aj 2 

15^bl 

+0.003 

Ay 2 

156bi 

+7.0 

At 

156bl 

+ 14.0 

2 At 

l60bl 

+0.65 

r l2 

l62bl 

+ 4.0 


l 64 bl 

+0.012 

n $ 

l66bl 

+10.0 


l6Sbl 

-.0 

/numbers 

17 Obi 

+0.5 

! 

I72bi 

+1.0 

$/hr (g d = 

174 bl 

+20.0 

k A# 

176b! 

+10.0 

m ( p g ) 

172>bl 

+10.0 


lSObl 

+3570.0 

t»lnl 

lS2bl 

+ 2397*0 

^2ml 

l 34 bl 

+2.0 

| & yiml 

l 36 bl 

+0.2 

|<f y 2m| 

lggbl 

+100 

m( Pay) 



+.0 


146 - 


Temporary storage 


190 bl 

192 bl 


TAPE NO. 2260 MO 

( 24 , 6 ) 

cl, c e 

t s 

c6, c s 

t e 
c s 
t a 

clO, a p 

c a 

t a 
s p 
1 n 

1 c 
i t 

8 p 

c2, s 1 

8 p 

c a 
t d 
s p 

8 p 


+.0 

F. 2 

CYCLE CONTROL PROGRAM 


x 2 0 
x 1 
o 

X 1 3 

x 2 

* 3 
a 1 

x 6 

* 3 

t 1 



} 

} 

: 


set run 
counter 
reset pass 
counter 

reset ordinate 
counter 
Read 1 Group 
save parts in 
11 b" registers 


a 156bl 
s 6 Obi 


c 2 

72 

d 1 



w 1 


Insert 2ftt 
in 60bl 

skip 1 Group fwd. 
read 1 Group fwd. 
save parts In 
"f" registers 

delay 
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o3, 

o4, 


o5, 

c6. 


C B 1 1 

S p 8 1 J 
8 p W 1 
8 p d 1 
c a x 1 3 | 
c p c 4 J 
s p k 1 
s p 1 c 4 

8 P P 1 

c a x 6 
td t 3 
8 p t 1 
C S 0 

- 

8 p 8 1 

c a x 1 6 
t d c 5 
s 1 72 
c a 0 
r c c 6 

a o c 5 

8 U X 1 4 

c p c 5 

8 p W 1 
a o x 3 
e p c 2 



skip back 
2 Groups 
delay- 

read 1 Group 
is pass counter 
positive? 
final comp. 


prelim, comp, 
save parts 
in "b" reg. 


skip back 
1 Group 


rerecord 
1 Group f-wd. 


delay- 

recycle if ord. 
counter is neg. 



8 i 4 0 S 


07. 


o9. 


ell, 


c a x 1 3 

c p c 9 

* 

a o x 1 

’ 

c p c 7 

8 p 2 1 

C 8 x 1 2 } 

s p a 1 

8 P C 2 

c a x 3 

8 U 1 

c p c 1 1 

a o x 1 3 

c s x 1 2 

8 p S 1 

C B X 2 

t 8 X 3 

a o x 3 



8 P C 1 0 
I N 

i c a 1 9 2 b 1 
i s p r 1 
i c a /2bl ) 


stop tape 

' 

is pass counter 
pos, ? 

9 

is run counter 

\ 

pos. ? 

final reprint 
skip back 
55 Groups 

, 

recycle 

* 

is ordinate 
counter + 21 


step pass counter 

( skip back 
53 Groups 
reset ordinate 
< counter to 

. < w> 

to 2nd pass 

delayed Print 

I-fS' 

. 

f delayed Print 
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c! 2 . 


cl3, 


i c a 72 bl 

1 s u 70 bl 
1 c p c 1 2 
i c a l 66 bl 
1 m r 170 bl 
its l 66 bl 

i s p c 1 3 

i a d I74bl 
i c p c 1 3 
i c a l 66 bl 
i d v 170 bl 
its l 66 bl 
i c a l 66 bl 
i s p r 1 
i c a 72 bl 

► 

its 70 bl 

i c a l 62 bl 
i t S i 91 b | , 
its 72 bl 

i c a 156 b!' 

► 

its 60 bl 


s p d 1 
c a x 5 
t d t 3 

s p t 1 

s p c 4 


kA$ 


c* 

k A $ + 30 

2<<= oc* 

delayed Print 
set up 
PXI/hr 
clear Z$/hr 
and H £ $ * 

Insert A t 

read 1 Group 
save parts 
in registers 


to Prelim. Comp 



F. 3 SUBROUTINE f 
DOLLARS PER HOUR SUMMATION 


1 t a f 2 
i c a 66 bl 
i a d 62 bl 
lad 26 bl 
its 52 bl 
i m r l64bl 
lad l 62 bl 
i m r 52 bl 
i a d 172 bI 
its 52 bl 
i a d 64bl 
lad 6 Sbl 
i s p 0 

F.4 SUBROUTINE g 
ELEVATION CHANGE LIMITER 


L - g x 

L - g x - g 2 = g d 
Sd 

+ g d 

m $ g a + n $ g d 2 
#/hr(g d =o)+m^g d +n $ g d 2 

$/hr(g d ) 

|/hr(g d )+ P ln 

#/hr(g d )+ P ln + P 2n =$/hr 


i t a g 5 
i c a 66 bl 
i c p g 3 
i s u 6 Sbl 
1 c p g 2 
i c a 6 Sbl 
1 s p g 5 
i c a 66 bl 


- 1 & VmJ 
I 
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1 

s p g 5 


e 3 , 

i 

a a 6 Sbl 

<ry + ISM 


i 

c P g ^ 



i 

c a 66bl 

SH 


1 

s P g 5 


g 1 *. 

i 

c s 6gbl 

S tjfrt 

85* 

i 

8 p 0 




F. 5 

SUBROUTINES h 



INCREMENTAL STORAGE, DAM 1 

hi, 

i 

t a h 3 



i 

a u 74 -bl 

y X - y 10 * y lB 


i 

t s 52bl 

y l8 


1 

m r gObl 

n ls y lB 


1 

a d 7^bl 

“lS + WlS 


i 

m r 52bl 

m ls y l8 + n lB y l0 2 


i 

a d 76bl 

a iO +m i 8 y iB +n i8 y iB 2 

h 3 * 

i 

s p 0 




INCREMENTAL STORAGE, DAM 2 

h2. 

i 

t a h 4 



i 

s u ll6bl 

y 2 *" y 20 “ y 2s 


i 

m r 120bl 

m 2s y 2s 


i 

a d llSbl 

s 20 + m 2a y 2s 

h 4 , 

i 

8 P 0 




F. 6 

SUBROUTINE k 



FINAL 

COMPUTATIONS 

kl» 

t 

a k g 



I 

N 
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1 c a 44bl f/>#/hr(t)/;y 1 (t)J = f64 

1 s u 54bl f64 - b64 

i d v 15 gbl /f64-b64/ / 2At 

Its 64bl = d P$/hr(t) 

3 dT T^Ttl 


1 s p k 2 
k2, 1 c a 46 bl 
i s u 56bl 
i d v 15Sbl 
Its 52bl 


space for occasional print out 
f/^$/hr(t-Hr 12 )/^y 1 (t)7 = f 65 
f65 - b 65 
ft6 5 - b 657 / 2At 

= d_ 3$/hr(t+T 12 ) 

at Txrti— 


1 8 p k 3 space for occasional print out 


k3, lea 52bl 

lad 64 bl 
1 s u 30bl 
1 s u 34 -bl 

1 d v 36bl 
1 n r i66bl 


p$/hr(t+r 12 ) 
at TT^Tt) 


d P $/hr(t) 

dt py^TtT 


+ d p$/hr(t+r 1 O 
57 ; y^t) 12 


-» P$/hr (t) + d__ P $/hr(t) + 


Ty 1 '('t) 


dt 3 ^TtT 


- P$/_hr(t) - p #/hr(t+T ) 

P y 1 (tT - ld 


P $/hr(t+r 12 ) 

-E 1 (t) / e-j^y-^t)) 


a yim 

= -E 1 (t) 


E-^t) 


_ oC 

i^Ty^TtTJ 

- o< _E^(t) = / y^(t) 


(y x (t)J: 


d_ p$/hr(t+T 12 ) 

dt a y-, (t) 


_ ir(t) + d_ 

dt 3 y 1 (t) dt 


1 d v 36b! 
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t s 66bl 

Sy 1 

ft) 

m r 66bl 

S7 1 

2 

m r 36bl 

8 1‘ 

ryi 2 

m r 36bl 

sS ± 

2 

a d I92bl 

1. 

z 


> y.is 

t s 192bl 

j 


c a lg4bl 

' 1 

t s 6gbl 


It 

c a b 1 

p T(IP 1 +i x ) 

a d 4-bl 

p rdP-L + i x ) v - 

c p k 4- 



s p g 1 

Limiter 

s p k 5 



c a 66bl 

s v 1 


a d 20bl 

s y x (4) + y 1 (t; 

t 8 2Obi 


It 

s p r 1 

delayed print 

c a 4-gbl 

fp$/hr(t)/3y ; 

8 U 5^1 

f67 - b67 

d v 15Sbl 

/f67-b677 / 2, 

a u 32bl 

- ] 

' 9 l/hr(t) 


i 

1 9 y 2 (t) 

d v 3Sbl 


i 


8 2 (y 2 (tJ) 

m r l66bl 

- 



I(IP 1 + i, f 


dt j ypTtl 


dt aypTtT 


4 

- 


= -E p ftJ 


Eg(t) 


Eg(t) 


2' J 2 


* Summation of all dam 1 penalty costs and spillage, for each 
variation described in 7.6a. 


** Similar sums for previous iteration. 
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i 

i 

i 

i 

1 

i 

i 

1 

1 

1 

1 

1 

i 

i 

k6, 1 

k7, i 
i 
i 
i 
1 
1 
1 

kg, s 


d v J>Sbl 


<=< _ 

s 2 (y 2 (t))£ 


t s 66bl / y 2 

m r 66 bl & y 2 2 

m r J&bl 

e 2 ‘y i 

m r 3Sbl J $ 2 2 

a d 192bl "] 

I< s 2 

t s 192 bl J 


c a lg6bl 

1 

t s 6Sbl 

ii 

c a 2bl 

p i (ip 

a d 6bl 

p i (ip. 

c p k 6 


s p g 1 

Limiter 


s p k 7 
c a 66bl 
a d 22bl 
t s 22bl 


SV 2 

y 2 (t) + s y 2 


E g (t) 


i 2 > 

+ 


8 p r 1 delayed print out 
c a 4b 1 
t s b 1 l 
c a 6bl 
t s 2bl j 
P 0 


= & y 2 (t) 


(ZP 2 + ij) 


set up violations counter 



-155 


F. 7 SIBROUTINE m 
PLANT 1 GENERATIONS 

ml, i t a m 6 

Its 50bl -Q 1 

i m r S6bl n it^l 

lad gLbl m.^ + 

i m r 5Obi ~ n lt Q]_ 2 

1 s u S2bl -y lt (^ 1 =°) - “ n it% 2 

lad 4 -gbl y 1 - y lt = h 1 

its 56bl h, 

i c a 4 gbl y^ 

i S U lCftbl y x - y lc = y 1L 

1 t s 52bl y 1L 

I c p m 2 

i m r lOSbl 

lad 106bl + n^y.^ 

1 m r 5 2bl miL y 1L + n^y^ 8 = J x 

1 s p m 3 

m2, 1 c a l6Sbl i x = 0 

m3, i t s 52bl 

lad 4 -bl save + £.P 
1 t s 4 bl 
i c a 52bl 

lad 50bl -Q x + i x = -Dj 
its 52bl 



i 

i 

i 

i 

i 

i 

i 

I 

i 

i 

i 

i 

i 

i 

i 

i 

i 

i 

1 

1 

1 

1 

i 

i 


m r 96 bl 

n ig D i 

a d 94bl 

“is + n ig D i 

m r 52 bl 

_m lg D l - n lg D l 2 

s u S 2 Sbl 

- h 10 - m lg D l - n lg D l 2 

a d 56 b! 

h l - h 10 - m Xg D l - n lg D l 2 

m r 92 bl 

T i^" h i “ h io “ m ig D i - n ig B i 2 -^ 

i r 52 bl 

/gi / d x 7 C-V-J = -g x 

t 8 62 bl 


c a 56 bl 

h i 

s u 9°bl 

h l " h ll = h ln 

t s 56 bl 

h ln 

m r 102 bl 

- n( Slm> h ln 

s u lOObl 

-m(g lm ) - n(g lB ) h ln 

m r 5 ^bl 

- n( BlB )h ln 2 

8 u 92>bl 

- m< Bl,n >h ln “ n( S 

t s 52 bl 


a d HObl 

g lm (ratca) - g^fhp 

c p m 4- 


c a 52 bl 


a p m 5 


c s HObl 

_g lm (ratea) 

s u 62 bl 

g l ~ g lm 

8 p U 1 1 

1 upper generation 

t s 64bl 1 

penalty cost 


g l/D 
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lea 62 b! ~g^ 

1 s p u 1 negative generation penalty cost 

lad 64 bl 

its 64 b 1 Ep x 

lea 17Sbl m(P, ) 

V 

its 50bl 
i c a 4 gbl j 1 
i s u lSObl y x - y lm 
i s p 2 1 

i a d 64 bl £ P x 

its 64 b 1 ' 

1 c a 7*H>1 mln 

1 8 u *1 min - *1 

i s p Z 1 

i a d 64 bl 

its 64 bl X 

i a d 4 bl KIPj +i x ) 

its 4 bl " 

m6, i s p o 

F. S SUBROUTINE n 
PLANT 2 GENERATIONS 

nl, 1 t a n 6 

its 5Obi -Qg 

i m r 126 bl n^Qg 

i a d 124b! 4 " n 2t^2 
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1 m r 50bl - m 2t^2 ~ n 2t^2 2 
i s u 122bl -Y2t^2 = °) ~ m 2t^2 ” n 2t^2 

lad 4 -gbl y 2 - y gt = h g 

1 t s 56 bl h 2 

i c a ^gbl y 2 

i a u iMbl y 2 - y 2o = y 2L 
i t a 5SD1 y 2L 
i c p n 2 

1 m r l 4 gbl n 2 L y 2 L 
lad l 46 bl m 2L + n 2L y 2L 
1 m r 52 bl m 2L y 2L + n 2L y 2L = i g 
1 s p n 3 

n2, 1 c 8 l68bl 0 

n3, 1 t s 52bl 

1 a d 6bl 51 ( 51 P 2 + J 2 ) 

1 t s 6bl " 

lea 52bl $2 
lad 50bl ~Q 2 + = —Dg 

Its 52bX -Dg 
1 m r 136bl +n 2g Dg 

lad 13 ^bl m 2g + n g g Dg 

1 m r 52bl - m 2g D 2 n 2g D 2 2 

i a u 128 bl -h 20 - m 2g D 2 - n^Dg 2 

1 a d 56b! h 2 - h 20 - m gg D 2 - n 2g D g 2 
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i m r I32bl 

r l 2 [h 2 - h 2Q - m 2g D 2 " n 2g D 2^ = g 2/D 2 

i m r 52bl 

Zg 2 / D 2^ /"- D 2^ = "g 2 

its 26bl 

-g 2 

led 56b! 

h 2 

i s u 130bl 

h 2 “ h 22 = h 2n 

its 56bl 

h 2n 

i m r l 42 bl 

- n(g 2 max ) h 2n 

i s u 14 -Obi 

~ m(g 2 max } - n(g 2m } h 2n 

i m r 5&bl 

- m( 8 2m > h 2n - n( > h 2 „ 2 

i s u 13Sbl 

-82m (h 2n =o) ~ ffi( «2m ,h 2n " n( «2m )h 2 n 2 

its 52bl 

-gg m (hg) 

i a d 150bl 

g 2m (ratea) - g 2m (hg) 

i c p n 4 


i c a 52bl 

-•WV 

i s p n 5 


n 4 , i c s 150bl 

-Bgmlrated) 

n5, i s u 26bl 

g 2 “ g 2m 

i s p u 1 j 


its 6Sbl 

I*p a 

i c a 26bl 

-g g 

i s p u 1 "] 


i a d 6gbl j 

’ Zp 2 

its 6Sbl j 

' 

i c a lggbl 

mtPgy) 



its 5Obi 
i c a M-Sbl y 2 
i s u lg 2 bl y 2 - y 2 max 
i s p Z I *\ 
i a d 6 Sbl 7 IP, 

its 6gbl J 

i c a ll6bl min 

i s u 4-gbl min - jfy ! 

i a p Z l 

i a d 6gbl IP, 
its 6$bl 

i a d 6bl £ (P 2 + i } ) 
its 6bl " 


n6, i a p 0 


pi, t a p 4 


F. 9 SUBROUTINE p 
PRELIMINARY COMPUTATIONS 


i c a l6Sbl 1 clear violations-count register 
1 t s 4 bl > 

its 6bl J 


i c a *!Obl fy 1 (t) 
i s u 50bl fy-^Ct) - by 1 (t) 

i d v 6Obi [fy 1 (t) - by^tj/At or 2*t = y x (t) 

its 44 b 1 y x {t) 



lea ^2bl 
i s u 52bl 
i d v 60 bl 

Its l4-6bl 

lea l6bl 
Its 66bl 
lea 20bl 
its 4gbl 
1 s p h 1 
Its 36 bl 
1 m r 44b 1 
its 5Sbl 
i s u S6bl 
i s p m 1 
i c a 62 b! 
its 32 bl 
i s p r 1 
i c a 64bl 
its 42bl 
i c a 22bl 
its 4gbl 

i s p h 2 
its Jgbl 


fy 2 (t) 

fy 2 (t) - by g (t) 

[fy 2 (t) - by g (t)]/4t or 2 At 

y 2 (t) 

L(t) 
y x (t) 
s 1 (y 1 (t)) 

8 1 (y 1 (t)) /Cl7 = 5 1 (t) 

S 1 (t) - F x (t) - -^.(t) 

-g 1 (t) 

save «.g 1 (t) 


delayed print out 
save Zp,(t) 


y 2 (t) 

s 2 (y 2 (t)) 
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i 

i 

i 

i 

i 

1 

1 

i 

i 

1 

1 

1 

1 

1 

1 

1 

i 

i 

i 

i 

1 

1 

i 


m 

r 

46 bl 

e 2 (y 2 (t) 

y 2 (t) 


t 

8 

54 b 1 

- S B <*> 



a 

d 

56bl 

S 2 (t) + 

S x (t) 


s 

u 

12bl 

S 2 (t) + 

$ 2 . (t) ■» F 2 (t) 


8 

P 

n 1 







gp(t) 



C 

a 

26bl 

d 



8 

P 

r 1 

delayed 

print out 


8 

P 

f 1 

$/hr(t) 



t 

8 

4 Obi 




8 

P 

r 1 

delayed 

print out 


C 

a 

40 bl 




a 

d 

72bl 

I t/hr(t) 


t 

8 

72bl 




c 

a 

20bl 

y x (t) 



a 

d 

112bl 

y 1 (t) + 



t 

8 

4 g>bl 

it 



8 

P 

h 1 

Si(yi(t) 

+ A!, i> 


m 

r 

44 bl 

8i(yi(t) 

+ A y x ) y ] _(t) = S 3 

L ( Vl +fly l ) 

t 

8 

34 b 1 

H 


8 

u 

Sbl 


+ A y^) - F 1 ( t ) = 

-Qi(t,yi + Ay 1 ) 

S 

P 

m 1 

ei ( Vi 

+ 4 7 !) 


8 

P 

f 1 

l/hrU^ + Ay 1 ) 


8 

u 

40 bl 

$/hr(t,3 

r l + ~ S/hr(t) 
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1 d V 112bl 

£^/hr(t,y 1 +Ay 1 ) - $/hr(t)]/ A y x 

its 30 bl 

= £$/hr(t )/& y 1 

i s p p 2 

space for occasional print out 

p2, lea 20bl 

yi<t) 

its 4gbl 

II 

lea 44bl 

yyt) 

lad ll4bl 

y x (t) + 

1 m r 3 6b 1 

8 l(y 1 (t)j * £ y x (t) + & yj = Si(^ 1 ^ 1 . + Ay i^ 

1 s u Sbl 

Sl (t lL y l + A y i) “ ^(t) = + A y x ) 

1 e p m 1 

g l (t ? y l + Ay l J 

1 8 p f 1 

$/hr (t j y^ + 4^) 

1 s u 40bl 

$/hr(t, y;L + A y^) - $/hr(t) 

1 d v ll4bl 

[$/hr(t, yi + A y^) - $/hr(t)]/A yi = J$/hr(t}/<* 

Its 24-bl 

ti 

lea 32 bl 

save -g 1 (t) 

Its 62 bl 


lea 4-2bl 

save ^P^Ct) 

its 64bl 


lea 22bl 

y 2 (t) 

lad 152 bl 

y 2 (t) + A y 2 

Its 4gbl 

H 

1 8 p h 2 

S 2 (y 2 + 

1 m r 46bl 

• • 

s 2 (y 2 +Ay 2 ) y 2 (t) = $ 2 (t,y 2 + Ay 2 ) 



i a d 56bl 

S 2 <V 2 *a y 2 ) + 

Si(t) 

i s u 12bl 

S 2 ( t ;y 2 + + 

5 x (t) - F 2 (t) = «Q 2 (t,y 2 

1 s p n 1 

g 2 (t,y 2 t« 2 ) 


i e p f 1 

|/hr(t,y 2 + ay 2 ) 


i s u 40 bl 

$/hr(t,y 2 +iy 2 ) - 

- $/hr(t) 

1 d v 152bl 

j$/hr(t,y 2 +Ay 2 ) . 

- $/hr(t)] / A y 2 

its 32bl 

= d |/hr/P y 2 


i c a 22bl 

y 2 ( t) 


its 4 gbl 



i c a 46 bl 

y 2 (t) 


i a d 154 -bl 

• • 

y2 + a y 2 


i m r 3Sbl 

e 2 (y 2 ) [ y 2 (t) + 

= $2 (t ’ y 2 + ^ y 2 ) 

i a d 5Sbl 

Si(t) + S 2 (Vy 2 

+ d’y 2 ) 

i s u 12bl 

s x (t) + 5 2 <t,*y 2 

+ A y g ) - F 2 (t) = -Q 2 (t,y 

i 8 p n 1 

+ A 


i s p f 1 

$/hr(t ? y 2 + A y g ) 


i s u 4 Obi 

^/hr(tjy 2 + A y 2 ) 

- $/hr(t) 

i d v 154 bl 

[$/hr(t,y 2 + d y 2 ) 

- $/br(t )]/4 y 2 

its 2 £>bl 

“ t$/hr/d y 2 


i c a 44 b1 

y x (t) 


i m r l60bl 

y 1 (t)f T i 2 ] 


i a d 20bl 

y x (t) + y 1 (t) [t 

121 * y l (t + T 12 ) 

its 4 Sbl 

y x (t + T i2^ 
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i c a 5Sbl 
i s u lObl 
i 8 p m 1 
i c a 46bl 
i m r l6obl 
lad 22bl 

its 4gbl 

i c a 54bl 
lad 56bl 
1 s u l4bl 
i e p n 1 
lea lgbl 
1 t 3 66bl 
1 s p f 1 
Its 40bl 
lea 54b1 
lad 34bl 
1 s u l4bl 

1 s p n 1 
1 8 p f 1 

1 s u 40bl 


S^(t) - §!(t + t 12 ) 

S x (t) - F x (t + t 12 ) = -(^(t + t 12 ) 

g x (t + T 1S ) 

y 2 <t) 

y a< t )[ T 12I 

y 2 (t)*y 2 (t)[r 12 ]-y g (t + t 12 ) 
y 2 (t + t 12 ) 

S 2 (t) 

$ x (t> + S 2 (t> 

Sj.Ct> + Sg(t) - F 2 (t + t 12 ) s -Q 2 (t + t 12 ) 
gg(t + t 12 ) 

L(t * t 12 ) 

|/hr(t + t 12 ) 

S 2 (t) 

$ 2 (t) + + 4 yi> 

S 2 (t) + 5 1 (t,y 1 + Ay : ) - Fg(t + t 12 ) 

5 -Q 2 (t + t 12 , y x + a y x ) 
g 2 (t + t 12 , y x n yi ) 

$/hr(t + t 12> y x + ay x ) 

$/hr(t + T12, yi + * yi ) - $/hr(t ♦ r 12 ) 
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i a v 112bl ^ $/hr(t + t 12 ) 

its 3^bl y x (t) 

space for occasional print out 
P 3 , i c a 44 bl y 1 (t) 

lad ll 4 -bl y^(t) + Ay^ 

1 m r 36b 1 s^^Cy^^Ct)) C y 1 (t) + A y^J * S 1 (t, y 1 + A y 1 ) 

lad 5^bl S x (t, y x + Ay x ) + S 2 (t) 

leu l 4 bl 5 x (t, y x +A^) + 5 2 <t) - Pg (t + r 12 ) 

* -^2^ + T 12* y l + Ay l J 
1 e p n 1 g 2 (t + t 12 ,■ y 1 + A y-j^) 

1 s p f 1 $/hr(t + t 12 , y x + 4 ^) 

1 8 u 4 -Obl #/hr(t + t 12 , y x + Ay 1 ) - $/hr(t + t 12 ) 

1 d v ll 4 -bl 2$/hr(t + t 12 ) 

1 t 8 26bl ^y-^t) 

p^, s p 0 

F. 10 SUBROUTINE d 
READ IN ONE GROUP 

dl, t a d 4 

c a d 5 1 set up first 

t d d 3 j address 

a 1 7 4 

d2, r d d 2 

d3, t s - - 


select tape 
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a o d 3 
s u d 6 
c p d 2 
s i k 0 g 
d 4 , s p o 
d5, t s b 1 
d6, t s 39bi 

sX| t fl 8 ^ 

t 8 X 1 1 

82 , 8 1 7 3 

a o x 1 1 

c p s 2 

8 l 4-og 

8 p W 1 
83, 8 p O 

wl, t a w 3 

C 8 X 7 

t s 2 7 
w2, a o 27 
c p w 2 


count 

stop tape 

let address 

last address 

P. 11 SUBROUTINE s 
SKIP BACK OVER N GROUPS 

set up counter 
pick up Bl. Marker 


stop tape 
delay 


F. 12 SUBROUTINE w 
DELAY 


W3, 8 p O 
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F. 13 SUBROUTINE t 
STORAGE TRANSFER 


tl, t a t 6 
c 8 x 1 & 
t S X 1 9 
c a x 4 - 
t d t 2 


I 


set up no. of 
registers transferred 
1 st Initial 
address 


t 2 , c a - - 
t 3 , t s - - 
t 5 , « 0 t 2 
a o t 3 
a o x 1 9 
c p t 2 
t 6 ,- s p 0 

F. 14 SUBROUTINE r 
DELAYED OUTPUT via MAG. TAPE 


rl. 


r 3» 


r 4 , 


i t 
0 U 

8 i 

S 1 

c a 
r c 
c a 
r c 
c a 


a r 2 
T 


' l 2 block markers 

7 0 J 

1 0 5 2^ 

r 3 


1 0 
r 4 
1 0 



5 0 


record 
the parts 
of m r a 


J 
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r5, r c r 5 

. stop in 

s i 7 0 

erased region 

b i 7 1 > 

s i 4 0 g J 

I N 

r 2 , i s p 0 

F. 15 SUBROUTINE z 
STORAGE ELEVATION PENALTY COST 

21 , i t a z 3 
1 c p 2 2 

its 52 b 1 violation y x 

i m r 52 bl y x 2 

i m r 50 bl m i y y x 2 or m 2 y y x 2 

i s p z 3 

z 2 , i c s l 6 gbl 0 
z3, i 8 p 0 

F. 16 SUBROUTINE u 
GENERATION PENALTY COST 

ul, i t a u 3 

i c p u 2 

its 52 bl violation g_ 

i m r 52 bl g x 2 

i m r 176 bl m g x 2 

i s p u 3 

u 2 , i c s l 6 Sbl 0 

u 3 , i s p 0 



F. 17 SUBROUTINE x 
SINGLE LENGTH PARAMETERS 


xl. 

+0 


run counter 

x2, 

* 5 

1 

"•■weeks - 1 

* 3 . 

+0 


ordinate counter 

x^, 

s i 

20 bl 

address of 1st register transferred 

*5, 

s 1 

40 bl 

fwd, new address 

x6, 

s i 

5 Obi 

bwd. new address 

x 7> 

s i 

4 0 

3*1 millisec. delay 

x 9» 

+ 1 

9 

"a words in group -1 

xio. 

+ 1 


no. block markers in 2 Groups -1 

xll. 

+0 


skip back counter 

xl2, 

+ 5 

2 

na.black markers in i + 1 Groups -1 

X13, 

-to 


pass counter 

xl 4 . 

c a 

39^1 

c a address of last word transferred 

xl6, 

c a 

b 1 

c a address of 1st word transferred 

xlS, 

+ 9 


# words transferred -1 

* 19 , 

+0 



x20, 

+ 9 


# runs «1 


START AT 

c 1 


F* l2 

CONTENTS OF TEMPORARY STORAGE REGISTERS 
IN THE ORDER OF INSERTION 


b 

1 

p KSp, 

+ i x ) 

previous step 

2 b 

1 

p £ (X P g 

+ ip) 

previous step 

4 b 

1 


+ i x ) 

current Step 

6 b 

1 

I<£p 2 

+ i 2 ) 

current step 


violation 

and 

spillage 

checks 
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gbl 

F x (t) 


1 0 b 1 

P l (t + T 12> 


1 2 b 1 

F 2 (t) 


1 4 b 1 

F 2^ + T 12^ 


1 6 b 1 

L(t) 


1 g b 1 

L(t + t 12 ) 


2 0 b 1 

y x (t) 


2 2 b 1 

*3 

CVI 


2 4 b 1 

C 6 4 

C 64 

2 6 b 1 

C 6 5, -C 4 o, ~c 4 i, -C 42 , -C43, «c 44 , 

-C45, C65 

2 g b 1 

C 6 7 

C67 

3 0 b 1 

C 6 2 

C62 

3 2 b 1 

C 6 6, -C 3 6 

G66 

3 4 b 1 


063 

3 6 b 1 

Sl (y 1 (t; )) 

s 1 < y“i ( t: )) 

3 g b 1 

e 2 (y 2 (t)) 

8 2 (y 2 (t)) 

4 0 b 1 

f y x (t), G 54 , G 57 


4 2 b 1 

f y 2 (t) , P-j^ ( t) 


4 4 b 1 

f 6 4 , Cl 


4 6 b 1 

f 6 5 , C 4 


4 g b 1 

f* 6 7 , y 


5 0 b 1 

b y x (t), ft ? m(P ly ) ; m( p 2 y^ 


5 2 b X 

b y 2 (t), y ls yiL,^x’" D i--ei m (h i ) -s 

ix^lr. e B ,0 



SrJgx , */ hr 


4 ' 
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5 4 b 1 b 6 4 , S 2 ( t) 

56 bl b 6 5 , h]_> h 2i 

Mbl b 6 7 , ^(t) 

6 0 b 1 At , 2At 

6 2 b 1 -C 3 $, -C 37 , -C3S, -C 39 

6 4 b 1 C6S 

6 6 b 1 L(t), £ y 1 , & y 2 

6 Sb 1 XP 2 

7 0 b 1 PZ#/hr 

7 2 b 1 £#/hr 
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Appendix G 

Examination of Round Off and Truncation Errors 
G.l Introduction 

In the application of the gradient method to the three 
plant problem, the simplest Lagrange -formula for approximating 
derivatives was used. This type of approximation is a source 
of inaccuracy that is fundamental to the method used, involv¬ 
ing both truncation errors because of the finite sampling and 
round-off errors because of the taking of differences of nearly 
equal quantities. On the other hand, the method operates to 
nullify rather than accumulate such errors in successive 
iterations, since each step evaluates anew the path of steep¬ 
est descent from the latest operating point. What is desired, 
then, is an examination of the order of magnitude of such 
errors in the trials made to establish that they are not 
excessive. 

G.2 Desired Accuracy of End Results 

The accuracy desired is dependent on the use to be made 
of the solution. In system operation, the inaccuracies in flow 
prediction will limit the desired accuracy of optimization. 

In project planning, uncertainties of other estimates will 
similarly make great accuracy unwarranted. In the following, 
rather severe requirements on the accuracy of the end result 
is assumed, and the consequent accuracy requirements on inter¬ 
mediate computations are examined. 
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In practice, the accuracy of each computation will 
fluctuate because of the variableness of the flow, load, and 
original drawdown curves. The objective of the following is, 
rather, to check order of magnitudes. Data obtained from the 
Whirlwind I digital computer is used to do this. 

In the case of dam I elevation ordinates are to be found 
to six decimal digits, ( 3552.12 ft.) so that with round-off, 
the elevations will be certain to the nearest tenth of a foot. 
The incremental storage of dam 2 is roughly ten times that of 
dam I. Hence the corresponding significant increment of 
elevation is one-hundredth of a foot. To allow for round-off, 
then, the elevations of dam 2 should be found to seven decimal 
digits ( 2 $$ 5 . 123 ). 

The changes in elevation to be made in the final steps 
of the iteration process should ordinarily be the same order 
of magnitude or less than the significant increments, i.e., 

0.1 ft. for dam I and 0.01 ft. for dam 2. Accordingly, both 
H and ‘Hi ,the changes in elevation at each step, should be 
found to two significant figures. 

G.3 Magnitudes of Terms in Euler Expressions 

The final equation, C71, given in Appendix E is: 
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The magnitudes of the terms in the bracket vary over a 
moderate range in normal operation" 1 " but vary over a much 
larger range when violations of operating limitations or 
spillage occur. From the data of the three plant problem, 
six sample weeks have been chosen to be representative. 
These sample weeks are chosen from the following periods: 

a) first ten weeks 

b) middle twenty-five weeks 

c) last sixteen weeks 

The values of the terms in (G-2) for these six weeks are 
tabulated in Table G-l. 

TABLE G-l 
COMPUTED VALUES 


1st : 

10 weeks 

middle 

25 weeks 

last 16 

weeks 

Weeks: a 

b 

c 

d 

e 

f 

9&/hir(t) ii inn 

Sift) 14 - 432 

46.3^5 

-0.591 

-0.634 

-0.744 

130.021 

?$/hr(i+ r a) negl. 
sy,(t) 

negl. 

negl. 

negl. 

negl. 

negl. 

A- iSj/hrCt) - 0.249 

23.439 

-1.133 

10.425 

-17.969 - 

-363.727 

jLyiikdtelW - 1.152 

d* diet) 

13.391 

-0.773 

5.923 

-15.377 - 

-171.661 


The magnitudes shown for the middle 25 weeks are representa¬ 
tive of normal operation during the drawdown season. The 
large positive values of 3 |/hr(t)/ jj t (t) during the first 
10 weeks are due to the fact that dams are near full in this 
interval and small explorations to higher elevations cause 

"^see figure 7.3 
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spillage and waste of water. (Thus, a safety margin is 
produced by decreasing elevation Ay, ft. below crest.) 

Large values are also found in the last 16 weeks period. 
However, these are due to severe violations of operating 
limitations (exceeding plant capacities) under the original 
operating mode. These large values reflect the necessity of 
relatively large corrections in this time interval to relieve 
the operating-limitation violations. 

G.4 Accuracy Requirements on Terms in the Euler Expression 
When taking the products of two numbers of differing 
accuracies, the more accurate number should be rounded off 
to contain one more significant figure than the least accurate 
factor. The total error is then taken to be the error of the 
least accurate factor. The final result is carried to the 
number of significant figures in the least accurate number 
(Scarborough, 43). Presuming << and $,(&(*)) to be avail¬ 
able to three significant figures the bracket 

/ ?t/hr (t) . 2$/hr(t+Tii) d$/kr(t) _ (t+fy) ) 

1 3Ut) M 3 i (t) J 

should be accurate to two significant figures » 

In the case of addition and subtraction of numbers of 
unequal accuracy, one should use in the numbers known to 
higher digits only one more decimal digit (tenths, hundredths, 
etc.) than is contained in the number known to the lowest digit 
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(Scarborough, 48 ). Presumably, some term in (G-2) may be 
known to only the minimum number of decimal places needed 
to obtain two significant figures for the bracket. The others 
then should have one additional decimal place. In line with 
this reasoning, we find the number of significant figures re¬ 
quired of each term in order that all have one more decimal 
digit (tenths, etc.) than the minimum. These are tabulated 
in Table G-II for each of the six weeks sampled for Tab]e G-I. 

TABLE G-II 

SIGNIFICANT FIGURES ASKED 

1st 10 weeks middle 25 weeks last 16 weeks 

Weeks: a b c d e f 

3 3 2 1 1 3 

1 3 3 3 3 3 

d Mr(Ur, x ) 2 3 2 2 3 3 

) y. Ci) 

Thus, if a standard is to be set, it would be desirable to 
obtain all terms of (G-2) to three significant figures (or at 
least have only few terms to only two significant figures.) 

G.5 The Time Derivatives 

Computations C6B and C69 are first derivatives of con¬ 
tinuous functions. As such, the computations involve both 
round-off errors and truncation errors. The latter results 
from the use of a "truncated” series to approximate the time 
derivative, instead of an infinite series. Both types of 
errors depend on the number of terms used in this approximation, 
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and on the spacing of the ordinates which sample the curve 
whose derivative is sought. 

To simplify programming and to save computer storage 

space, the Lagrange type of formula involving ordinates was 
chosen rather than the Newton type of formula involving dif¬ 
ferences. Since it can be shown that the Newton and Lagrange 
interpolation formulas are basically identical, it follows 
that the truncation error made by using either will be the 
same if both employ the same number of terms. To further 
simplify programming, expressions using equally spaced 
ordinates were chosen. 

An examination of differentiation formulas for equally 
spaced ordinates reveals advantage for those having an odd 
number of points and the derivative obtained at the midpoint. 
Not only is the formula simpler than others of the same degree, 
but also the truncation error is smaller. Such symmetrical 
formulas are thus to be preferred for all but the end points 
of the time interval being optimized. (Milne, 40, pp 96-99). 
The derivatives with respect to time have, except at the end 
points, been approximated by the simple expression: 

h = 'zhl'**' ^ (0-3) 

where h is the spacing between ordinates */<-, ty, Hz, and the 
derivative is being evaluated at the y, ordinate.^" This ex¬ 
pression is the three point Lagrange formula for central 

Itfis here considered to be a general variable, not necessarily 
storage elevation. 
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derivatives. The remainder term or truncation error for this 
expression can be shown to be (40, p. 96 ): 



where jj is the value of the third derivative of f/ for some 
value of the independent variable between ij 0 and y 2 . 

G.6 Estimates of Truncation Errors 

As is evident from (G-4), the truncation error is depend¬ 
ent on the higher derivatives or high frequency content of the 
function being differentiated. This will naturally vary with 
different stream flow, load, and original drawdown curves. 
Moreover, it will vary markedly at different steps of the 
solution by the gradient method. As shown in figures 7.5b and 
7.6, the gradient method successively smoothes the cost curves, 
as the solution progresses. In fact, the solution tends 
towards a flat cost curve, as far as resource and operating 
limitations will allow. Therefore, the truncation errors de¬ 
crease as the solution is approached, which is of course very 
desirable. 

The data thus far obtained from Whirlwind gives successive 
reductions in cost without reaching a point of diminishing 
returns or otherwise indicating that a minimum has been ap¬ 
proached. Therefore, the initial data will be examined, with 

relatively large-magnitude high-frequency components in the 
cost curves, to obtain order of magnitudes of truncation 
errors at the beginning of the smoothing process. 
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The cost function is known only at the sampled points. 
Hence an exact determination of the third derivative in 
(G-4) cannot be obtained. A reasonable approximation can be 
obtained, however, under the assumption that the derivative 
sought does not change rapidly in the short time interval 
between ij Q and . Under such conditions, an approxima¬ 
tion to the (ntl)th derivative is given by the formula (Milne, 

40 , p.128): 


(»*') (»+i) /n+i\ rn* n (H*t) , 

h n (z) = -(. )% + (< ij^.t + 


(G-5) 


where ("f ) , (”2 0 > etc. are the binomial coefficients 

belonging to the exponent n*l; Z is a value of the independent 
variable between the ordinates y o and . This value of Z 
is not necessarily equal to the value of the independent var¬ 
iable occuring in the error formula (G-4), but is useful under 
the above mentioned assumption that the derivative sought 
changes slowly in the interval considered. 

The best data on hand gives the values of the first deriv¬ 
atives 

jsL 3f>/kr Ci) j__ H/kr ft* Till (G _ 6) 

elt 9 i, (t) Mt 2 y, (t) 

Therefore, equation (G-5) is used to obtain an estimate of the 
third derivative as the second derivative of the computed first 
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derivatives. Using (G-5): 

h‘ f (n cz) = y, - (?) h, + (I) 

* % - 1 y. + s. 


Therefore: 


f“\z) 


hi - +i 


(G-7) 


(G-S) 


M/hrCt+fc) 

The values of £ y i Ct) and 9 $,(t) are 

relatively smooth in each of the three regions selected above, 
but change more abruptly at the boundaries of these regions. 

To better examine the range of the corresponding third deriva¬ 
tives, they have been computed by equation (G-3) and tabulated 
in Table G-III for sampled weeks near the beginning and the 
center of each region.^ 


TABLE G-III 
THIRD DERIVATIVES 

1st 10 weeks middle 25 weeks last 16 weeks 


2$/h \r (t) 
dt 3 By.Ct) 


beginning 

0.0002 


jdl 0.0758 

di* Bid) 


center 

+.0646 

0.3315 


beginning 

0.1950 

0.0225 


center 

0.2372 

0.1151 


beginning center 
0.0290 0.3274 

0.0455 1.556 


In Table G-IV are shown the corresponding values of C63 and C 69 
plus or minus their probable truncation error with one week 


^Extreme variations in C63 and C69 were selected and were not 
necessarily chosen from the same weeks. 
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spacing of ordinates. Using the same values for the third 
derivatives, the corresponding truncation errors with h equal 
to one-half week are shown in Table G-V. 

TABLE G-IV 


TIME DERIVATIVES ± ESTIMATED TRUNCATION ERRORS 

WEEKLY SPACING 


1 st 

10 

beginning 

, 

dt 

17.11 

9 t/kr Ct) 

*i Ct) 

± 0.0016 

dt 

15. 83 

Mr(ur it ) 

act) 

± 0.619 

-L. \y 

weeks 

center 

5.901 

t 0.5275 

3.717 

j- 

0.312 

middle 

beginning 

6.677 

- 1.592 

1.152 

+ 

0.133 

25 

weeks 

center 

9.31$ 

f 1.937 

1.370 


0.939 

last 

16 

weeks 

beginning 

17.96 

i 0.236 

1.587 


0.371 

center 

31.61 

t 6.757 

103.1 

4- 

12.707 


TABLE G-V 

TIME DERIVATIVES± ESTIMATED TRUNCATION ERRORS 
HALF WEEK SPACING 

* 

a Mm 

dt Jy.ct) 


1 st 

beginning 

17.11 

4 

0.0004 

15.83 

•+ 

0.15 

10 

weeks 

center 

5.901 


0.13 

3.717 

± 

0.078 

middle beginning 
oc 

6.677 

+■ 

0.39 

1.152 

♦ 

0.045 

weeks 

center 

9.31$ 

+ 

0.4$ 

1.370 

♦ 

0.235 

last 

beginning 

17.96 

4 

0.059 

1.587 


0.093 

16 . 

weeks 

center 

31.61 

4- 

1.68 

103.1 

+ 

3.17 
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TABLE G-VI 

ESTIMATES OF NUMBER OF SIGNIFICANT FIGURES 
OBTAINABLE IN EARLY ITERATIONS 
WITH ESTIMATED TRUNCATION ERRORS 

C68 C 69 

h * 7 h= 3.5 h * 7 h= 3.5 

1st beginning 5523 

10 

weeks center 1222 

middle beginning 1223 

25 

weeks center 1212 

last beginning 3322 

16 

weeks center 2223 

Considering truncation errors alone, and assuming the 
truncation errors of tables G-IV and G-V to be precise, the 
number of figures to which C68 and C 69 are known, with an 

accuracy of better than one-half the highest digit, are tab¬ 
ulated in table G-VI, for one week and half-week ordinate 

spacings. It is seen that with one week spacings the origi¬ 
nal desire that all terms in the bracket be known to two or 
three significant figures is not met at the start of the 
smoothing process, whereas with half-week spacing this desire 
is met. Remembering that this evaluation is for the beginning 
of the smoothing process, and expecting the higher derivatives 
to be reduced at each step of the iteration, it is concluded 
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that the order of magnitude of error is not excessive and that 
one weeks spacing for this data is probably appropriate as far 

as truncation error is concerned. 

G.7 Higher Order Lagrange Formulae 

To examine the variation of truncation error as the order 
of the Lagrange formula increases, consider the next higher 
symmetric differentiation formula, the five point formula: 

. . r "l , M (s) 

tf. +j fr * ( 0 - 9 ) 

To evaluate this error, we obtain from (G-5): 

h'f C¥, (z) *&-(?]& + (’)y, * 

G-10) 

= & - H - H + if. 

Hence 

h y cz)~y 4 -4y s t -4* + y. (0 -ii) 


The fifth derivatives have been computed using (G-ll) at an 
ordinate in each of the three regions considered above. The 
corresponding truncation errors are tabulated in table G-VII. 
Comparison of these errors with those of table G-V indicates 
that the truncation errors in the five point formula are, in 
these cases, larger than those in the three point formula, 
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when one week spacing is used. Thus, with one week spacing, 
five points are spread over too large an interval, with ex¬ 
cessive variations in the derivatives (C68, C 69 ) therein. 


TABLE G-VII 

Weeks Spacing 

Truncation Errors with 5 Point Formula 



A 3$/hr Ct\ 

At W.tt) 

A, 9 $/Wfr-fc) 

At WXt) 

1st 10 wks (center) 

.3778 

.2057 

Middle 25 wks (center) 

-.9484 

-.4715 

Last 16 wks (center) 

-11.685 

-7.979 


6.8 Round Off Error 

To examine roundoff errors in obtaining C68 and C69, 
information is needed on the magnitudes of ZSl h di L 

/ » i.ct) 

and &t/hr(t+Tn.)JA V>(1) . Values of these quantities were 

obtained in computer runs made on June 16, 1952. Plant char¬ 
acteristics, load, and flows used in this run were the same 
as those in the runs considered above. However, the dams 
were operated with full reservoirs instead of being drawn- 
down. Nevertheless, the data obtained makes possible an 
estimate of roundoff error. In table G-VIII are tabulated 
values of $/hr(t), - 

M.ct) ; y,(t) aMt) Wit*) 

l.tArLft + Kil and —i ilkJAi r nl - for every fifth 

W.(t) i*tt) 

week taken from the data of June 16. This provides a good 
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comparison of effects of each change on the $/hr cost of 
operation. Note that all are positive because 

under this operation any increase in elevation increases 
spillage. The last three weeks sampled are during the 
spillage season. Proposed operation would have overloaded 
machines, and the derivatives are governed almost entirely 
by penalty functions which dictate a reduction in loading. 

In all cases it is seen that the term 

3 $/hr(t + Vi) 

is negligible in comparison to some other terms which enter 
into the bracket, eq. (G-2), in spite of the fact that eleva¬ 
tions were at crest height, where variations in elevation are 
more influential due to the imminence of spillage. 

Again, values are selected from the three periods used 
above, sampling data at weeks 2, 6, 16, 2$, 36 and 45. These 
are tabulated in table G-IX along with the number of signifi¬ 
cant figures needed in each in order that the bracket, eq, 
(G-2) be accurate to two significant figures. 

In all the time derivatives, we have used the simple form 

1 0j- (G-12) 

2 At 

If the resultant is to be known to q significant figures, the 
numerator should be known to q figures and the denominator to 

q*l, Considering the denominator to be exact, we require only 
that the numerator be known to q significant figures. 
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TABLE G-VIII 

COMPUTED VALUES AT SAMPLE WEEKS 



$/hr (t ) 

3t/hr(t) 

nfhr(t) 

2$/hr(t) 

Week 


a h, ct) 

a J, ct; 

J HxCt) 

5 

1332.983 

58.654 

3855.13 

42.423 

10 

2138.885 

61.222 

4014.68 

45.068 

15 

2048.016 

60.995 

4009.32 

44.018 

20 

2162.194 

62.705 

4123.39 

44.987 

25 

2506.848 

65.801 

4313.75 

48.828 

30 

2343.302 

63.880 

4189.05 

46.834 

35 

2078.47 

61.449 

4039.56 

44.263 

40 

43,282.18 

11.974 

889.66 

-15,415.03 

45 

.567 x io 7 

-.247 X 10 6 - 

-.265 X 10 S 

2.0833 

50 

.127 X 10 7 

.48177 

76.432 

.223 X 10 1 


n/kM 

9$/lr(t+fu) 

P $>/kr(t+T/z) 

Week 

2 % (t) 

3 y, ct) 

9 if, ct; 

5 

9453.98 

-.6 X 

-24 

10 

1780.529 

10 

9927.41 

-.1 X 

10” 23 

1863.91 

15 

9792.48 

-.6 X 

10" 24 

1333.95 

20 

10,035.73 

-.1 X 

10-23 

1888.35 

25 

10,643.96 

-.1 X 

10" 23 

1999.61 

30 

10,343.34 

-.1 X 

io“ 23 

1936.56 

35 

9361.16 

-.6 X 

10-24 

1840.52 

40 

-.45 X 10 7 

.0358 

-.252 X 10 7 

45 

1364.53 

.6250 

275.520 

50 

.342 X 10 g 

.0651 


.634 X 10 7 


Data of June 16, 1952 
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TABLE G-IX 


VALUES AND SIGNIFICANT FIGURES 
ASKED OF TERMS IN EULER EXPRESSION 



)$/kr(t) 
3 y,c t) 


A 

dt 3 vxt) 


Week 

Sig. 

Value Figures 

Sig. 

Value Figures 

Value 

Sig. 

Figures 

2 

53.7711 

3 

10.26 3 

21.7$ 

3 

6 

59.4771 

3 

6.32 2 

3.79 

2 

16 

61.0901 

3 

1.35 2 

0.549 

1 

23 

65.6234 

3 

-3.19 2 

-1.46 

2 

36 

60.3227 

3 

-7.47 2 

-3.12 

2 

45 

247,020.6 

3 

-110,300 3 

1.225 

0 


Data of June 16, 1952 


The differences corresponding to the numerator of (G-12) 
are tabulated in table G-X for the selected weeks. Using the 
required number of significant figures, tabulated in table 

G-IX for xn 4 J_ jSArfor*) 

dt ) <!, (t) H 3 if, (t) 

one then obtains the number of significant figures needed in 

1 t/ktt) w a 

3}, Or) 3 jf, Ct) 

These are also tabulated in table G-X. To illustrate the 
procedure used, consider week two. From table G-X, in week 2, 

d 2&Mt) 

dt ? (t) 

is needed to three significant figures. It follows that the 








TABLE G-X 


COMPUTED VALUES AND ACCURACIES ASKED OF PARTIAL DERIVATIVES 
TO LIMIT ROUNDOFF ERROR IN TIME DERIVATIVES 


Week 

2 

6 

16 

23 

36 

45 

$/hr{ t) 

1305.37 

1972.245 

2040.173 

2497.963 

1956.699 

5,670,040 

9Wlnr(i) 

3 i, C*) 

3576.603 

3901.254 

4017.607 

4301.489 

4003.333 

-26,516,305 

A 2t/tr(i. )_ J&/h 

3 i c-t) ncutu) 3 ict-it) 

413.756 

S3.290 

13.901 

-44.702 

-104.744 

-1,551,157 

§ sig, fig. needed in d$/kr(£). 

3 i,(+) 

4 

5 

4 

5 

4 

5 

4 

5 

3 

4 

4 

5 

d$/hr(t+17x) 

9 i(t) 

1540.194 

1312.195 

1337.329 

1992.692 

1793.393 

275.521 

$t/kr(t+V>+U) Btlhr(t+r a -dt) 

304.513 

53.03 

7.70 

-20.427 

-113.847 

17.133 


4 

4 

4 

4 

3 

0 

9 it*) 

5 

5 

5 

5 

4 

1 
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difference: 

a a/hdtl = sjMtif) - l 

(0-13) 

is also needed to three significant figures. Examination of 
the relative magnitudes of 


A 3 t/kr(t) 

3 t 


and 


s id) 


in table G-X shows that the less accurate term in (G-13) 
should be known to four figures and the other should be known 
to five figures. 

It is recalled that these requirements stem from the 
rather severe requirement that elevations be accurate to 
better than the nearest tenth and nearest hundredth 
at dams 1 and 2 respectively. Lesser requirements on the 
accuracy of the end results will naturally need less accuracy 
in the intermediate computations. 


G.9 Other Errors 

Truncation errors involved in the difference-method of 
evaluating 

tJt/kmsL and iMdhfl 

3 Ct) 9 4cct) 

are not treated, because such errors are eliminated by the 
methods of 6.2 and 6.3 using analytic expressions for these, 
partials. 

Further evaluations of the accuracy requirements of 
plant characteristics are being made but are not ready for 
this publication. Indications are, however, that this is not 
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a serious problem unless very accurate results are required. 
Rather, the problem is to determine how much simpler the 
characterization can be made and still obtain a desired accuracy. 

G.10 Conclusions 

As stated previously, such error analysis can only provide 
estimates of errors and guides for future work. The examination 
presented indicates that truncation and round off errors in the 
application considered are not excessive, and no fundamental 
limitation that would prevent the use of the method has appear¬ 
ed. Adjustment of sample spacing and modification of the degree 
of the Lagrange formulas for different problems should of 
course be investigated. In the use of the method on an actual 
system, only periodic checking of errors will insure satis¬ 
factory operation. 
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COMPUTER PROGRAMS 

Program. A program is a sequence of actions by which a computer 
handles a problem. The process of determining the sequence of actions 
is known as programming. 

Flow diagrams. A flow diagram is a series of statements of what the 
computer has to do at various stages in a program. Lines of flow indi¬ 
cate how the computer passes from one stage of the program to another. 

Coded program . Programs and flow diagrams are largely independent 
of computer characteristics, but instructions for a computer must be ex¬ 
pressed in terms of a code. A set of instructions that will enable a com¬ 
puter to execute a program is called a coded program, and the process 
of preparing a coded program is known as coding. 

Orders and operations . Individual coded instructions are known as 
orders and call for specific operations such as multiply, add, shift, etc. 

The computer code. The computer code described here is that of 
Whirlwind I, an experimental computer using binary digits, single¬ 
address order code, parallel operation, and electrostatic storage. It is ex¬ 
pected that computers of this type will ultimately achieve an average 
speed of 50,000 operations per second. 

COMPUTER COMPONENTS 

Registers and words . A register has 16 digit positions each able to 
store a one or a zero. A word is a set of 16 digits that may be stored in a 
register. A word can represent an order or a number. 

Arithmetic element . Arithmetic operations take place in the arith¬ 
metic element, whose main components are three flip-flop registers, the 
A-register, the accumulator, and the B-register (AR, AC, BR). The 16 
digit positions of AR starting from the left are denoted by ARO, ARl, . . 
AR15. Similarly for AC, BR. Words enter AC through AR; BR is an 
extension of AC. 

Storage . The term "register" by itself refers to the main electrostatic 
storage, which consists of 2 11 or 2048 registers, each of which is identi¬ 
fied by an address. These addresses are 11-digit binary numbers from 
0 to 2047. The computer identifies a register by its address. 

Input-output. All information entering or leaving the computer is tern - 
porarily stored In the input-output register (IOR). The computer regu¬ 
lates the flow of information between the internal storage and IOR, and 
also calls for any necessary manipulation of external units. The descrip¬ 
tive names of the input-output orders were chosen for photographic film 
reader-recorder units, but the orders are applicable to other types of ex¬ 
ternal equipment. 

Control element . The control element controls the sequence of com¬ 
puter operations and their execution. Instructions are obtained from stor - 
age in the form of individual orders, each of which is represented by a 
single word. 

Inter-connections . The four main elements (storage, control, arith¬ 
metic, and input-output) are connected by a parallel communications sys¬ 
tem, known as the bus. 

REPRESENTATION OF ORDERS 

Operation section . When a word is used to represent an order the first 
(left-hand) 5 digits, or operation section, specify a particular operation 
in accordance with the order code. 

Address section . The remaining 11 digits, or address section, are 
interpreted as a number with the binary point at the right-hand end. In the 
majority of orders this number is the address of the register whose con¬ 
tents will be used in the operation. In orders si, sr, the number specifies 
the extent of a shift; in rf, rb, the number specifies an external unit; in ri, 
rs, the address section is not used. 

Example . The order ca x has the effect of clearing AC (making all the 
digits zero) and then putting into AC the word that is in the register whose 
address is x. If q is a quantity in some register, the order needed to put q 
in AC is not caq but cax, where x is the address of the register that con¬ 
tains q. 


REPRESENTATION OF NUMBERS 

Single-word representations . When a word is used to represent a 
number the first digit indicates the sign and the remaining 15 are numeri¬ 
cal digits. For a positive number the sign digit is zero, and the 15 numeri¬ 
cal digits with a binary point at their left specify the magnitude of the num - 
ber. The negative -y of a positive number y is represented by comple¬ 
menting all the digits, including the sign digit, that would represent y. (The 
complement is formed by replacing every zero by a one and every one by 
a zero.) In this way a word can represent any multiple of from 
2“^ - 1 to 1 - . Neither +1 nor -1 can be represented by a single 

word. Zero has two representations, either 16 zeros or 16 ones, which 
are called +0 and -0 respectively. 

Overflow — increase of range and accuracy . With single-word 
representation the range is limited to numbers between - l and 

I - 2“ 1 ^ . Programs must be so planned that arithmetic operations will not 
cause an overflow beyond this range. The range may be extended by using 
a scale factor, which must be separately stored. Accuracy can be in¬ 
creased by using two words to represent a 30-digit number. 

COMPUTER PROCEDURE 

Sequence of operations . After the execution of an order the program 
counter in the control element holds the address of the register from which 
the next order is to be taken. Control calls for this order and carries out 
the specified operation. If the order is not sp or cp(-) the address in the 
program counter then increases by one so that the next order is taken from 
the next consecutive register. Thesp and cp(-) orders permit a change 
in this sequential procedure. 

Transfers . A transfer of a digit from one digit position to another af - 
fects only the latter digit position, whose previous content is lost. 

Negative zero . The subtraction of equal numbers produces a negative 
zero in AC, except when AC contains +0, and -0 is subtracted from it. 

Manipulation of orders. Words representing orders may be handled in 
the arithmetic element as numbers. 

Procedure in the arithmetic element . The execution of an addition in¬ 
cludes the process of adding in carries; this process treats all 16 digits 
as if they were numerical digits, a carry from AC 0 being added into AC 15. 
A subtraction is executed by adding the complement. Multiplication, divi¬ 
sion, shifting and round-off are all executed with positive numbers, com¬ 
plementing being performed before and after the process when necessary. 
For round-off the digit in BR 0 is added into AC 15. 

NOTATION FOR CODING 

Addresses . A coded program requires certain registers to be used for 
specified purposes. The addresses of these registers must be chosen be¬ 
fore the program can be put into a computer, but for study purposes this 
final choice is unnecessary, and the addresses can be indicated by a sys¬ 
tem of symbols or index numbers. 

Writing a coded program . Registers from which control obtains or¬ 
ders may be called action registers, and should be listed separately from 
registers containing other information, which may be called data regis¬ 
ters. A coded program is written out in two columns; the first contains 
the index number of each action or data register, and the second column 
indicates the word that is initially stored in that register. In many cases 
part or all of a word may be immaterial because the contents of the regis¬ 
ter in question will be changed during the course of the program. This 
state of affairs is indicated by two dashes, for example, ca —. 

The abbreviations RC, OR . Abbreviations used in referring to the 
register that contains a certain word or to the word in a certain register 
are 

RC . . . = (Address of) Register Containing . . . 

CR . . . = Contents of Register (whose address is) . . . 

The symbol ri x . When an address forms part of an order it is repre¬ 
sented by the last 11 digits of a word whose fir st 5 digits specify an opera¬ 
tion. An address x that is not part of an order is represented by the last 

II digits of a word whose first 5 digits are zero, which is equivalent to 
specifying the operation ri. Thus the word for an unattached address x 
may be written ri x. It could also be written x x 2“ 1 ^ . 



THE ORDER CODE 


AC = Accumulator 

x is the address of a storage register; 


AR - A-Register BR - B-Register 
n is a positive integer; k designates an external unit 



Operation 



Order 

Name 

Code 

Function 


Decimal 

Binary 


ri — 

read initially 

0 

00000 

Take words from external unit until internal storage is full. 

rs -- 

remote unit stop 

1 

00001 

Stop external unit. 

rf k 

run forward 

2 

00010 

Prepare to use external unit k in forward direction. 

rb k 

run backward 

3 

00011 

Prepare to use external unit k in backward direction. 

rd x 

read 

4 

00100 

Transfer to register x a word supplied by external unit. 

rc x 

record 

5 

00101 

Arrange for transfer of contents of register x to external unit. 

ts X 

transfer to storage 

8 

01000 

Transfer contents of AC to register x. 

td x 

transfer digits 

9 

01001 

Transfer last 11 digits from AC to last 11 digit positions of register x. 

ta x 

transfer address 

10 

01010 

Transfer last 11 digits from AR to last 11 digit positions of register x. 

cp(-)x 

conditional program 

14 

onio 

If number in AC is negative, proceed as in sp; if number is positive disregard the cp(-) 
order, but clear the AR. 

sp X 

subprogram 

15 

01111 

Take next order from register x. If the sp order was at address y, store y + 1 in last 11 
digit positions of AR. 

ca x 

clear and add 

16 

10000 

Clear AC and BR, then put contents of register x into AC. If necessary, add in carry from 
previous sa addition. 

cs X 

clear and subtract 

17 

10001 

Clear AC and BR, then put complement of contents of register x into AC. If necessary, 
add in carry from previous sa addition. 

ad x 

add 

18 

10010 

Add contents of register x to contents of AC, storing result in AC. 

SU X 

subtract 

19 

10011 

Subtract contents of register x from contents of AC, storing result in AC. 

cm x 

clear and add magnitude 

20 

10100 

Clear AC and BR, then put positive magnitude of contents of register x into AC. If neces¬ 
sary add in carry from previous sa addition. 

sa x 

special add 

21 

10101 

Add contents of register x to contents of AC, storing result in AC and retaining any over¬ 
flow for next ca, cs, or cm order. Only orders 1 through 15 may be used between the sa 
order and ca, cs, or cm orders for which the sa is a preparation. 

ao x 

add one 

22 

10110 

Add the number 1 x 2 _1 5 to the contents of register x. Store result in AC and in register x. 

mr x 

multiply and round off 

24 

11000 

Multiply contents of register x by contents of AC; round off result to 15 numerical digits 
and store in AC. Clear BR. 

mh x 

multiply and hold 

25 

11001 

Multiply contents of register x by contents of AC and retain the full product in AC and the 
first 15 digit positions of BR, the last digit position of BR being cleared. 

dv x 

divide 

26 

11010 

Divide contents of AC by contents of register x, leaving 16 numerical digits of the quotient 
in BR and to in AC according to sign of the quotient. (The order si 15 following the dv order 
will round off the quotient to 15 numerical digits and store it in AC.) 

si n 

shift left 

27 

11011 

Multiply the number represented by the contents of AC and BR by 2 n . Round off the result to 
15 numerical digits and store it in AC. Disregard overflow caused by the multiplication, but 
not that caused by round-off. Clear BR. 

sr n 

shift right 

28 

11100 

Multiply the number represented by the contents of AC and BR by 2~ n . Round off the result 
to 15 numerical digits and store it in AC. Clear BR. 

sf X 

scale factor 

29 

11101 

Multiply the number represented by the contents of AC and BR by 2 sufficiently often to make 
the positive magnitude of the product equal to or greater than 1/2. Leave the final product in 
AC and BR. Store the number of multiplications as last 11 digits of register x, the first 

5 digits being undisturbed. 


NOTES ON THE ORDER CODE 

Effect of operations . The functions of the various orders are de¬ 
scribed above. It is to be assumed that AR, AC, BR, and the register 
whose address is x are undisturbed unless the contrary is stated. 

AR. AR is primarily a buffer register for passing words into AC. 
After orders cax, csx, adx, sux, sax, and aox it contains the number 
originally contained in register x. After orders cmx, mrx, mhx, and 
dv x it contains the magnitude of the contents of x. The effect of sp x and 
cp{-)x is stated above. No other order changes the contents of AH. 

BR . A number stored in BR always appears as a positive magnitude, 
the sign of the number being assumed to be that indicated by the sign 
digit in AC. This convention has no effect on the logical result of the 
operations involving BR except that when BR contains a number that will 
be used later it is necessary to retain the appropriate sign digit. 

Alarms . If the result of an arithmetic operation exceeds the register 


capacity (i.e., if overflow occurs), a suitable alarm is given except as 
mentioned in connection with orders sax and sin. 

Shift orders . A multiplication overflow in si is lost without giving 
an alarm, but an overflow from round-off gives an alarm. Orders srO 
and slO only cause round-off, an alarm being given if an overflow occurs. 
The integer n is treated modulo 32, i.e., sl32 = slO, sl33 = sll, etc. 

Scale factors . If all the digits in BR are zero and AC contains ±0, the 
order sfx leaves AC and BR undisturbed and stores the number 33 in 
the last 11 digit positions of register x. 

Division . Let u and v be the numbers in AC and register x when the 
order dvx is used. If ju j < |v| the correct quotient is obtained and no 
overflow can arise. If |u|>|v| overflow occurs and gives an alarm. If 
u = v/ 0 the dv order leaves 16 ones in BR and round-off in a subsequent 
si 15 would cause overflow and give an alarm. If u = v = 0 a zero quotient 
is obtained. 







REVISIONS IN THE WHIRLWIND I ORDER CODE 


Ootober, 19^9 through May, 1951 



Operation 


Order 

Name 

Dec. 

Code 

Binary 

Function 

A. Permanent addition to the 

Code: 



ck x 

check 

11 

01011 

Stop the computer and ring an alarm if contents of register x is not Identical 
with contents of AC; otherwise proceed to next order. 

B. Minor changes not affecting previous 
functions of certain orders: 

sa x 

special add 

21 

10101 

Add contents of register x to contents of AC, storing result in AC and retaining 
any overflow for the next C£l, £s, or cm. Only orders 1 through 15, 23, 30, 
and 31 may be used between sa and the ca, cs, or cm for which the sa is a pre¬ 
paration. Use of any other operation 'Between sa and ca, cs, or cm will result 
in the overflow being lost completely, with no other eTfe'cT on t'Re normal func¬ 
tion of the intervening operation. 

sl*n 

shift left without 
roundoff 

27 

11011 

Multiply contents of AC and BR by 2 n , as in si n, but do not roundoff nor 
clear BR. 

ar*n 

shift right without 
roundoff 

28 

11100 

Multiply contents of AC and BR by 2* n , aB in sr n, but do not roundoff nor 
clear BR. 

C. Temporary order likely to become 

permanent: 

qe x 

exchange 

13 


Exchange contents of AC with contents of register x (original contents of AC 
to register x, original contents of register x to AC). 

D. Major changes in the order 

code: 



rs— 

(remote unit) stop 

X 

00001 

To be discontinued: at present it has the temporary function of stopping the 
computer. 

rf k 

run forward 

2 

00010 

Discontinued. Replaced by operation si. (described below). 

rb k 

run backward 

3 

00011 

Discontinued. Replaced by operation si (described below). 


E. Tentative new order (not yet completely 
defined): 


si k 

select in-out unit 



Select the particular piece of terminal equipment (with mode and direction of 





operation, if necessary) specified by the address. 


P. Present temporary orders (to be discontinued 
when replaced by the more general ajL, rd 
and rc operations): 


qh x 

h-ttxia eet 

6 

00110 

Transfer contents of AC to register x; set the horizontal position of all 
aisplajr scope beams to correspond to the numerical value of the contents 
of AC. 

qd x 

D-scopes display 

7 

00111 

T * a ?!! fe ^ co ? tentB of AC to re 8ister xj set the vertical position of the beams 
of the display scopes to correspond to the numerical value of the contents of 

AC; display (by intensifying) a spot on the face of the D-diBplay scopes. 

qf 1 

P-scope display 

23 

10111 

Same as operation £d, except display a spot on the face of the P-display 
scopes. 

qr r> 

read/ shift right 

30 

11110 

Perform two logically distinct functions: l) Cause the photoelectric tape 
input reader to read the next character cont/ilnlng a hole in position 7 from 
tape into digits 0 through 5 of FF Register #3 (holes or no holes in tape po¬ 
sitions 1 to 6 becoming ones or zeros in digit columns 0 to 5 respectively). 

2) Shift the contents of AC and BR to the right n times. The sign digit is 
shifted like any other digit and zeros are introduced into the left end. 

(no roundoff, no BR clear, no sign control). Note: If more than 3 milli¬ 
seconds (« about 60 orders) elapBes between one qr and the next, the reader 
stops. This must not happen except where 3" of blank tape has been provided 
for the purpose. 

qr # n 




Same as qr n above, except the mechanical tape reader, rather than the photo¬ 
electric tape reader, is caused to read. The mechanical reader can read line- 
by-line (i.e. no blank tape is required at places where the reader is allowed 
to stop). 

qp n 

punch/ shift right 

31 

11111 

Perform two logically distinct functions: 1) Cause the paper tape output 
equipment (punch or printer or both, depending on the settings of the 
switches) to record ope character with holes or no holes in positions 1 to 

6 on tape corresponding respectively to ones or zeros in digit columns 10 
through 15 of AC, and with a hole In position 7. 2) Shift right as in op¬ 

eration jjr. 

qp*n 




Same as qp n above, except no hole is put in position 7. 

qs n 

index camera 

___ 

12 

01100 

Perform two logically distinct functions: 1) Move the next frame of film into 
place in the display scope camera and open the shutter If it is not already 
oDen. The shutter is closed manually when the film is removed for developing. 

2) Shift right as in operation jr. 


■•The fact that the six largest binary digits of the address section of orders si n, sr n , qr n 
and qp n are normally zero (and in any case are disregarded by the step counter which counts 
the shifts) is exploited to permit the addition of an extra variant, as described, to each of 
the orders mentioned. The star (*) in sl*n , sr*n, qr*n , and qp*n implies that digit 6 (the 
2 y - 512 digit of the address) will be a one, while in ~si n , etoT, with no star, digit b is 
to be a zero. This can be accomplished during preparation of tape by typing si (800 + n) 
for #l*n. 


























